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ABSTRACT 


Coherent  speckle  noise  is  modeled  as  a  multiplicative  noise  process  that  has  a 
negative  exponential  probability  density  function.  Using  a  homomorphic  transfor¬ 
mation,  this  speckle  noise  is  converted  to  a  signal-independent,  additive  process. 
The  speckled  images  are  randomly  jittered  from  frame-to-frame  against  a  uniform 
background  to  simulate  image  motion  and/or  platform  jitter. 

Multiple  images  are  logarithmically  transformed  and  ensemble  averaged  in  the 
bispectral  domain.  The  bispectrum  ignores  this  image  motion  so  no  blurring  results 
from  the  ensemble  averaging.  Object  Fourier  magnitude  and  phase  information 
are  also  retained  in  the  bispectrum  so  that  the  resultant  image  can  be  uniquely 
reconstructed.  This  value  is  then  exponentiated  to  complete  the  image  reconstruc¬ 
tion  process. 

Since  speckle  masks  the  resolution  of  details  in  the  noisy  image  and  effectively 
destroys  the  object  structure  within  the  image,  it  is  seen  that  image  reconstruction 
using  bispectrum  estimation  results  in  images  that  regain  their  object  structure. 

Both  one-dimensional  and  two-dimensional  images  were  tested  using  separate 
bispectral  signal  reconstruction  algorithms  for  each. 
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1.0  INTRODUCTION 


When  an  object  is  illuminated  by  a  coherent  source  of  electromagnetic  radiation, 
and  this  object  has  a  surface  structure  that  is  rough  on  the  order  of  a  wavelength  of 
the  incident  radiation,  a  speckle  pattern  results.  A  fully  developed  speckle  pattern 
will  appear  chaotic  and  disorganized.  Even  when  the  speckled  object  is  imaged,  the 
presence  of  speckle  noise  is  seen  as  a  collection  of  spots  superimposed  on  the 
actual  image.  Thus,  speckle  is  considered  as  noise  that  degrades  an  image  whereby 
information  present  in  the  image  is  masked.  This  noise  can  be  present  in  either  the 
pupil  plane  (farfield)  or  the  image  plane. 

Studies  related  to  the  occurrence  of  speckle  phenomenon  are  not  new.  Dainty 
(19S4)  relates  a  brief  history  of  observations  of  speckle  patterns,  beginning  with 
Newton's  attempt  to  explain  star  twinkling.  Exner.  in  the  late  19th  century, 
sketched  a  speckle  pattern  created  by  a  candle  seen  through  fogged  glass.  Dehaas. 
in  the  early  20th  century,  actually  photographed  such  a  pattern.  In  the  early  1960's, 
with  the  advent  of  the  laser,  the  investigation  of  this  phenomenon  accelerated. 
Coherent  speckle  phenomenon  can  occur  is  such  diverse  imagery  as  synthetic  aper¬ 
ture  radar  (SARi.  acoustic  imagery,  ultrasound,  x-ray  scattering,  electron  scatter¬ 
ing.  microwaves,  and  obviously  laser  illumination. 

What  really  is  important  is  how  speckle  degrades  an  image.  Speckle  increases  the 
size  of  the  minimum  resolution  patch  obtained  with  a  given  aperture  as  compared 
to  the  same  size  aperture  using  incoherent  illumination  (Kozma  and  Christensen. 
1976).  The  resolution  of  the  image  is  reduced  making  certain  feature  identification 
difficult.  Tine  Maximum  Likelihood  Estimate  (MLE)  of  a  speckle  image  is  obtained 
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from  ensemble  averaging.  This  technique  is  a  highly  effective  way  of  removing 
speckle  noise  (Sadjadi.  1990;  Jain  and  Christensen,  1980;  Lim  and  Nawab,  1981). 
Multiple  images  are  necessary  to  perform  this  procedure.  Perfect  registration  of 
each  image  with  the  other  is  required  so  that  the  resultant,  averaged  image  is  not 
blurred.  This  concern  often  precludes  using  ensemble  averaging  due  to  the  diffi¬ 
culty  in  registering  images  exactly.  This  is  especially  true  when  the  image  is  chang¬ 
ing  its  position  from  frame-to-frame. 

A  signal  processing  technique  called  bispectrum  estimation  can  perform  ensemble 
averaging  on  multiple,  shifted,  speckled  images.  The  bispectrum  is  defined  in  this 
paper  to  be  the  Fourier  transform  of  the  triple  correlation  of  an  image.  Image 
reconstruction  with  the  bispectrum  has  the  following  beneficial  properties: 

1.  The  bispectrum  is  insensitive  to  linear  phase  shifts; 

2.  The  bispectrum  retains  both  Fourier  magnitude  and  phase  information: 

3.  The  bispectrum  is  insensitive  to  certain  additive  noise  processes. 

This  technique  is  a  spectral  domain  technique  as  opposed  to  a  spatial  domain 
technique,  where  much  of  the  current  digital  speckle  reduction  techniques  operate. 
This  thesis  will  review  the  statistics  of  speckle,  speckle  as  noise,  and  how  this  noise 
affects  the  image.  A  thorough  examination  of  speckle  computer  modeling  will  be 
conducted  and  several  methods  of  generating  speckled  objects  digitally  will  be  re¬ 
viewed  in  detail  with  examples  of  speckle  imagery  given.  I  will  also  review  the 
bispectrum  and  how  it>  use  with  homomorphic  processing  will  aid  in  image  recon- 


struction  of  speckle-degraded  images.  Two  separate  bispectral  reconstruction  algo¬ 
rithms  are  used  for  one-dimensional  and  two-dimensional  images. 


2.0  BACKGROUND 


The  technique  I  will  use  to  reconstruct  the  image  of  a  speckled  object  will  be  a 
wholly  digital  technique.  By  this  I  mean  that  each  speckled  image  will  be  in  a 
digitized  format  and  will  be  operated  on  totally  using  a  computer  on  a  pixel-by¬ 
pixel  basis.  I  state  this  at  the  beginning  so  as  to  differentiate  this  technique  from 
physical  methods  of  reducing  speckle.  For  example,  McKechnie  (1984)  lists  sev¬ 
eral  methods  one  can  do  to  reduce  speckle  before  it  is  imaged,  such  as: 

1.  Illuminating  the  object  with  a  temporally  partially  coherent  source: 

2.  Illuminating  the  object  with  a  spatially  partially  coherent  source; 

3.  Observ  ing  the  speckled  object  through  a  moving  aperture  that  performs  time¬ 
overaging  as  it  moves  across  the  pattern: 

4.  Observing  the  pattern  through  a  finite  aperture  that  acts  as  a  low-pass  filter. 

Guenther,  et  a!.  (1978)  has  one  of  the  earliest  papers  listing  several  different  digi¬ 
tal  filters  that  could  be  applied  to  reduce  speckle.  They  describe  the  basic  ensem¬ 
ble  averaging  technique  (where  several  images  are  averaged  together  to  reconstruct 
the  image)  as  well  as  the  spatial  processing  technique  using  an  averaging  window 
(which  moves  across  the  image  spatially).  Two  simple,  digital,  non-linear  filters 
(square-root  and  squaring)  are  also  described.  Jain  and  Christensen  (1980)  review- 
similar  techniques  as  well  as  a  homomorphic  Wiener  filter  for  digital  speckle  re¬ 
moval.  Lim  and  Nawab  (1981)  compare  other  digital  techniques  for  speckle  reduc¬ 
tion  such  as  low-pass  filtering  in  the  frequency  (Fourier)  and  density  (logarithmic) 
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domains  as  well  as  a  method  they  call  the  “short  space  spectral  subtraction  image 
restoration  technique.”  They  discuss  briefly  the  homomorphic  approach  to  image 
restoration  using  a  filter  built  specifically  for  this  method  but  do  not  describe  this 
filter  in  any  detail,  though  they  do  discuss  the  benefits  of  using  the  homomorphic 
approach.  Sadjadi  (1990)  reviewed  the  before-mentioned  averaging  techniques,  a 
median  filter,  local  statistical  filters,  an  adaptive  filter,  and  a  sigma  filter  and 
compared  these  with  the  homomorphic  approach  to  speckle  reduction.  The 
homomorphic  approach  allows  additive-noise  reduction  techniques  to  be  applied  to 
multiplicative  noise  conditions.  A  better  description  of  a  local  statistics  filter  used 
in  homomorphic  processing  is  given  by  Arsenault  and  Levesque  (1984).  A  general 
description  of  signal-dependent  noise  is  presented  in  this  article,  also. 

A  comparison  of  phase-retrieval  algorithms  was  performed  by  Fienup  (19S2)  and 
some  aspects  of  these  algorithms  are  used  to  restore  speckled  images  (Idell,  et  al. 
19S7).  Cederquist,  et  al,  (1988)  used  a  phase  retrieval  algorithm  for  farfield,  com¬ 
puter-generated  speckle.  These  algorithms  are  useful  because  phase  information  is 
lost  when  spectral  processing  a  recorded  image  using  the  power  spectral  density 
(PSD).  Recovering  the  phase  helps  one  reconstruct  the  image  uniquely.  Bispectra! 
processing  retain>  the  phase  of  an  object  and  it  is  this  idea  that  makes  it  a  possible 
technique  to  perform  image  reconstruction. 

Kuan,  et  al.  (1987).  derived  an  adaptive  restoration  filter  for  speckled  images.  This 
article  also  brings  up  some  important  points.  One  key  point  is  the  difference  be¬ 
tween  a  noisy  object  and  a  noisy  imaging  system.  A  noisy  object  occurs  when  an 
incident  coherent  wavefront  is  scattered  by  the  object  such  that  random  interfer- 


ences  occur  among  the  dephased,  but  still  coherent,  reflected  waves.  This  is  the 
speckle  case  I  am  investigating.  The  noisy  imaging  system  results  from  a  randomly 
variable  transmitting  medium,  such  as  the  atmosphere,  and  this  is  the  case  for 
stellar  speckle.  Instead  of  being  object-related,  it  is  medium-related,  i.e.,  a  point  is 
individually  speckled  throughout  the  image  due  to  system  effects.  Crimmins  (1985) 
used  a  geometric  filter  on  speckled  images  and  then  compared  his  results  with  a 
3X3  median  filter  used  on  the  same  images.  Safa  and  Flouzat  (1988)  used  mor¬ 
phological  technique?  on  remotely-sensed  speckled  SAR  images. 

All  of  these  previous  noise-reduction  methods  operate  on  the  speckled  image  in 
the  spatial  domain  (except  for  the  homomorphic  Wiener  filter).  Recently, 
Marathay,  et  al.  (1989)  used  computer-simulated  speckle  patterns  and  applied 
third  and  fourth  order  intensity  correlations  to  restore  the  image.  The  speckle  pat¬ 
terns  generated  were  not  imaged  but  rather  found  in  the  pupil  plane.  Newman  and 
Van  Yracken  (19S9)  combine  the  bispectrum  with  photon-bias  correction  tech¬ 
niques  (see  Northcott.  et  al,  1988)  to  recover  multiple,  shifted  objects  in  a  uniform 
background  of  photon  noise  (e.g.,  Poisson).  The  signal-to-noise  ratio  (SNR)  of 
these  objects  and  background  noise  was  less  than  one!  Freeman,  et  al.  (198$).  used 
the  bispectrum  on  one-dimensional  (1-D),  infrared,  astronomical  speckle  data  to 
recover  the  phase  of  a  stellar  object.  Ayers,  et  al,  (198S)  compare  the  triple  corre¬ 
lation  with  the  Knox-Thompson  phase  retrieval  technique  for  astronomical  speckle 
and  their  result  indicate  that  the  triple  correlation  is  more  robust  for  recovering 
phase  than  the  Knox-Thompson.  Again,  as  a  reminder,  astronomical  speckle  is 
different  than  the  speckle  considered  in  this  thesis.  Astronomical  speckle  can  also 


be  termed  “incoherent  speckle”  while  the  speckle  modeled  in  this  thesis  is  coher¬ 
ent.  Mavroidis,  et  al,  (1990)  found  that  the  bispectrum  could  not  be  used  for  recov¬ 
ering  coherently  imaged  objects  through  atmospheric  turburlence.  This  was  based 
on  the  theory  of  phase  closure  and  its  relation  to  the  bispectrum.  What  they  did  not 
do  was  perform  the  homomorphic  transformation  before  they  attempted  image 
recovery.  As  will  be  shown  later,  this  is  a  key  parameter  to  successfully  reconstruct 
an  image. 

The  goal  of  this  thesis  is  to  use  recently  developed  algorithms  for  bispectrum  esti¬ 
mation  to  reconstruct  multiple,  shifted  images  degraded  by  speckle  noise.  This  is  a 
conceptual  stud}  to  determine  if  the  bispectrum  is  robust  enough  to  handle  the 
image  degradation  resulting  from  speckle.  The  aim  is  to  implement  solutions  from 
an  engineer's  viewpoint  so  that  practical  solutions  can  be  obtained.  Therefore,  in 
some  cases,  a  median  filter  may  be  applied,  after  bispectra]  averaging  and  the 
image  reconstruction  has  been  implemented,  to  further  smooth  the  image. 


3.0  SPECKLE 


To  properly  model  speckled  images,  I  must  attempt  to  review  the  underlying  phys¬ 
ics  and  mathematics  behind  speckle  theory.  Since  speckle  is  random,  it  lends  itself 
quite  well  to  statistical  analysis  and  the  literature  is  very  complete  regarding  this 
aspect  of  speckle.  By  using  key  concepts  from  this  area,  it  will  be  shown  that 
speckle  obeys  certain  constraints  that  are  conducive  to  computer  modeling. 

3.1  Electromagnetic  Field 

The  general  method  I  will  use  to  describe  the  electric  field  vector  of  the  propagat¬ 
ing  source  is  defined  by  several  authors  (Gaskill,  Haus)  and  is  reviewed  here. 
Obviously  any  representation  of  the  electric  field  vector  must  also  obey  Maxwell’s 
equations,  which  in  turn  leads  to  the  wave  equation: 

v:E  -  pc  =  0  (1) 

3  t  ■ 

where: 

E  is  the  electric  field  intensity  vector  E(x,y,z;t) 

(i  is  the  permeability 

e  is  the  dielectric  constant 

,  3:  32  3: 

v  -  is  the  Laplacian  operator:  - -  + - 7  +  - 7  • 

3  x :  3  y *  3  z- 

One  solution  to  the  wave  equation  is  found  by  writing  the  vector  E  as  (with 
r  =  (x.y.z)  being  a  vector  component): 

E  [ r :  t  ]  =  At  r  i  e.\p[  j  2tt  vt ] 
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where: 


A(r)  is  the  amplitude  function 
v  is  the  mean  frequency. 

In  imaging  I  am  concerned  primarily  with  spatial  quantities  and  not  temporal  ones. 
Thus,  I  can  suppress  the  time  dependence  in  the  previous  equation.  Furthermore, 
since  the  electric  field  is,  in  general,  a  complex  quantity,  I  can  express  it  as: 

A  (r )  =  Re j.Af  r  i]  +  jlm  fA  (r)]  =  |  A  ( r  )|  exp  [j<J>(r)]  (3) 


where: 

Re  denotes  the  real  part  of  the  field 
1m  denotes  the  imaginary  part  of  the  field 

IAU)  I  denotes  the  magnitude  of  the  field 

cj>(  r )  denotes  a  phase  function. 

Therefore,  this  is  the  representation  of  the  complex  amplitude  of  the  wave  field.  I 
will  use  this  notation  throughout  the  thesis.  Another  quantity  of  interest  will  be  the 
intensity,  defined  as: 

1(E)  =  |  A  (  r  )  |:  (4) 

It  is  the  squared  modulus  of  the  complex  field.  Implicit  in  the  definition  of  the 
complex  field  that  is  a  solution  to  the  wave  equation  is  that  it  is  assumed  to  be 
perfectly  monochromatic  and  linearly  polarized.  Linear  polarization  is  required  so 
that  each  component  can  be  treated  as  a  separate  random  walk  (Goodman.  19S6) 
and  so  that  each  component  will  be  statistically  independent.  Monochromaticity  L 
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required  because  if  other  wavelengths  are  present  then  each  one  can  produce  its 
own  speckle  pattern  causing  the  overall  statistics  of  the  speckled  image  to  be  pre¬ 
cluded  by  the  effect  of  all  the  individual  statistical  distributions  for  each  wave¬ 
length.  The  reflected  field  remains  linearly  polarized  and  monochromatic,  though 
it  does  undergo  a  uniformly  distributed  phase  change  upon  reflection. 

3.2  Object  Surface  Statistics 

Once  this  electromagnetic  field  strikes  the  object,  it  is  reflected  from  this  object. 
The  surface  is  assumed  to  reflect  all  of  the  energy  and  any  loss  is  due  to  the 
destructive  interference  of  the  reflected  waves.  It  is  at  this  point  that  the  speckle 
phenomenon  occurs.  The  derivation  for  the  object  field  statistics  that  I  will  use  is 
the  classic  random  walk  (Beckmann  &  Spizzichino,  1964;  Goodman,  1976,  19S4, 
1985,  19$6j.  This  approach  requires  some  assumptions  as  follows  below. 

A  “rough  surface’'  is  a  surface  which  scatters  the  energy  of  the  incident  plane 
wave  into  various  and  random  directions  upon  reflection.  The  “smoothness”  of  this 
surface  is  dependent  upon  the  relative  wavelength-to-surface-roughness  relation¬ 
ship.  If  this  relationship  is  of  the  order  of  a  wavelength  (where  a  2r r  phase  excur¬ 
sion  can  occur  that  is  uniformly  distributed),  speckle  can  occur.  If  the  phase  is  not 
uniformly  distributed  then,  even  with  the  presence  of  a  large  number  of  scatterers 
(a  central  limit  requirement  [Stark  &  Woods, 1986]),  the  speckle  intensity  will  not 
have  its  required  negative  exponential  distribution.  Similarly,  a  large  number  of 
scatterers  are  required  so  that  Gaussian  random  variables  can  be  used  as  the  field 
amplitude  model.  Shadowing  effects  (where  a  portion  of  the  surface  “shades"  an 
adjacent  scatterer)  and  multiple  scattering  from  the  surface  are  neglected:  only  one 
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reflection  of  the  incident  field  occurs  at  the  surface  before  it  interferes  with  an 
adjacent  coherent,  reflected  wave.  The  density  of  the  scatterers  (number/m2)  is  not 
considered  and  there  is  always  enough  return  signal  to  measure  adequately,  i.e, 
unconstrained  intensity  is  assumed  so  as  to  preclude  Poisson  statistics 
(photocounts).  I  assume  that  this  surface  reflected  field  (or  its  analytic  representa¬ 
tion)  can  be  modeled  as  a  circular,  complex,  Gaussian  random  variable  whose 
probability  density  function  is  given  as: 


p(Ar,  A, 


exp 


. 


~(  ( Af  )2  +  (  A;  )2  ) 

2tto2 


(5) 


where: 

Aris  the  real  part  of  the  field,  A(r) 

Aj  is  the  imaginary  part  of  the  field,  A(r) 
o  is  the  standard  deviation. 


This  is  a  complex  Gaussian  process  where  both  components  are  independent,  iden¬ 
tically  distributed  with  zero  means  and  equal  variances.  By  using  a  probability 
transformation  one  can  show  that  (Goodman,  19S4)  the  modulus  of  this  Gaussian 
field  is  Raxleigh  distributed: 


P  (  A(  r ) 


(6) 


where: 

| At  r)  |  =  \  (  Ar  )'  +  ( A j ) ‘ 
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It  can  be  shown  that  the  intensity  1  =  |  A(r )  J“  (Goodman,  1984)  has  a  negative 
exponential  probability  density  function  given  by: 


(7) 


From  the  previous  development,  2o2  =<|ak|:  >  ,  where  <  >  denotes  the  ensemble 

average.  Notice  that  the  value  within  the  <  >  is  simply  the  individual  phasor  intensi¬ 
ties.  therefore  1  can  rewrite  Equation  7  in  the  more  familiar  form: 


Pd  ) 


1 

<]> 


(8) 


Another  parameter  is  called  the  contrast  of  the  speckle  pattern  and  is  defined  as  an 
inverse  SNR.  For  fully-developed  speckle,  the  contrast  should  equal  one  (Good¬ 
man.  19S5:  Cederquist,  19SS).  The  contrast  is  given  by: 


All  of  these  first-order  statistics  will  be  important  when  I  create  computer¬ 
generated  speckle  because  they  will  be  used  to  verify  that  the  speckle  patterns  are 
actually  “speckled." 

3.3  Propagation  And  Imaging  Of  The  Speckle  Field 

Since  I  am  dealing  with  images,  their  representation  must  of  necessity  be  real  and 
non-negative  because  image  intensity  (or  its  field  modulus)  can  never  be  less  than 
zero.  Second-order  statistics  will  help  in  describing  the  image  formation  process  in 
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this  case  because  these  statistics  define  the  spatial  structure  of  the  field  (Dainty, 
1976).  I  have  reviewed  several  cases  and  will  relate  key  results  here.  Goldfischer 
(1965)  derived  the  autocorrelation  (ATC)  and  the  PSD  for  a  uniform  diffuse  object 
viewed  in  the  farfield  at  the  pupil  plane,  thus  not  imaged.  The  farfield  is  that 

ttL2 

region  a  distance  z  from  the  object  that  satifisfies  the  constraints  of  z  »  - 


(Gaskill),  where  L  is  the  size  of  the  aperture  and  X  is  the  wavelength  of  the  inci¬ 
dent  radiation. 


The  farfield.  or  Fraunhofer,  diffraction  pattern  intensity  is  proportional  to  the  com¬ 
plex  ATC  of  the  object  and  is  thus  a  measure  of  the  spatial  frequency  spectrum  of 
the  power.  In  other  words,  the  ATC  and  the  PSD  are  Fourier  transform  pairs 
(Yariv).  Goldfischer's  model  used  a  uniform,  and  hence  a  statistically  stationary 
object.  Since  images  are.  in  general,  non-stationary,  it  will  be  shown  that  his  re¬ 
sults  are  a  special  case  of  the  more  usual  case  of  non-uniformity.  Enloe  (196") 
later  extended  these  results  to  an  imaged,  uniform  object. 

There  exists  a  Fourier  transform  relationship  between  the  speckle  field  and  its 
farfield  propagation,  given  by  (Enloe): 

CC 

Aj  ( £.  T) )  =  / ’/a n  (x,y)  exp  \ (x  t  +  yp  )1  dx  dy  (10) 

where: 

A0(x,y)  is  the  object  speckle  field 
Aj(  c.p  )  js  the  farfield  distribution 


At  the  farfield  an  aperture  is  located  at  the  pupil  plane.  This  aperture  could  be  a 
pupil  stop  some  distance  from  a  lens  or  the  lens  itself.  The  lens  then  focuses  this 
coherent  radiation  onto  a  detector  located  at  the  image  plane.  The  resultant  field  at 
the  image  plane  is  given  by: 

00  r  -n  1 

A2(v,u>)  =  //H(&,Ti)A1(£,Ti)exp  U^-(v£+wti)  d£  dq  (11) 

-00  L  -1 

w here: 

H(Cr|)  is  the  System  Transfer  Function  (Fourier  transform  of  pupil  func¬ 
tion) 

f  is  the  focal  length  of  the  lens 

v.w  merely  denote  a  different  coordinate  system  from  that  of  the  object. 

The  lens  acts  as  a  Fourier  transformer  (Goodman,  1968).  Thus  the  field  at  the 
image  plane  is  the  Fourier  transform  of  the  field  in  front  of  the  lens,  modified  by 
the  aperture.  The  intensity  at  the  image  plane  is  found  from  Equation  11  times  its 
complex  conjugate  (denoted  by  “)  by; 

I  ( v.  to  )  =  A:  (  v.  w  )( A;  ( v,  to  ) )  *  (12) 

Notice  that  the  intensity  at  the  image  plane  is  nothing  more  than  the  autocorrela¬ 
tion  of  the  object  convolved  with  the  imaging  system. 

3.4  Nonstationary.  Second-Order  Image  Statistics 

The  previous  development  was  for  stationary  objects.  Lowenthal  and  Arsenault 
(1970)  studied  the  case  of  an  imaged,  nonstationary  object  being  speckled.  Their 
development  is  important  because  several  important  relationships  regarding  the 
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object  structure  and  the  second  order  statistics  of  the  speckle  field  arise.  To  begin 
with,  the  object  f(r)  is  composed  of  two  components:  t(r)  is  the  actual  object  func¬ 
tion  that  acts  to  modulate  the  rough  surface,  which  is  denoted  by  d(r).  This  non- 
uniform  object  is  represented  by  (neglecting  the  vector  representation  for  r  for 
notational  simplicity): 

f(r)  =  t(r)d(r)  (13) 

The  image  plane  field,  g(r),  is  related  to  the  object  field  via  the  convolution  opera¬ 
tor  (*)  with  the  system  impulse  response,  h(r): 

g(r)  =  f(r)*h(r)  (14) 

Relating  d(r)  to  the  statistics  discussed  previously,  d(r)  is  stationary  and  a  circular, 
complex  Gaussian  random  variable  with  independent  real  and  imaginary  parts  hav¬ 
ing  zero  means  and  equal  variances.  The  function  can  also  be  thought  of  as 
wideband,  random  noise  (Lim  and  Nawab,  1981).  An  important  point  here  is  that, 
even  though  the  object  is  non-stationary,  its  noise  function  is  stationary  and  thus, 
for  all  cases  of  Ur),  we  can  find  the  parameters  of  g(r).  Referring  to  Equationl4.  1 
can  express  the  mean  of  the  image  as  (h(r)  is  the  same  function  throughout  the 
ensemble): 

<  g(r)  >  =  <  f(r)  >  *  h(r)  (15) 

Since  <  d(r)  >  =  0.  then  <  f(r)  >  =  <  t(r)  >  <  d(r)  >  =  0.  and  therefore  <  g(r)  >  =  0. 
The  variance  of  the  amplitude(or  field)  is  given  for  real,  non-negative  images  as 

'■s  ^  — 

(from  the  standard  relationship  o'  =  <x‘>-<x>'  (Friedan): 

(ot, (r))'  =  <  |  g ( r)  | *  >  -  |<g(r)>|*  (16) 
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Since  the  image  mean  is  zero  and  the  magnitude  squared  of  the  amplitude  equals 
the  intensity,  I  can  rewrite  the  variance  of  the  image  amplitude  as: 

(ag(r))2  =  <  I  (r)  >  07) 

This  interesting  result  says  that  the  variance  of  the  image  amplitude  is  equal  to  the 
mean  of  the  image  intensity.  For  the  variance  of  the  intensity  ]  can  write: 

(Oj  (r)) 2  =<I  (r):>  _  <  I  (r)  >2  08) 

The  definition  of  a  statistical  ATC  function  is  given  by  Goodman(19S5)  for  the 
intensity  (since  intensity  is  real)  as: 

Ri i ( r i .  r 2  '  =  <  I  (  r  j  )  I ( r 2  )  >  09) 

This  function  allows  us  to  compare  the  intensities  at  r  j  and  r2  over  the  ensemble. 
By  letting  r  {  =  r->  =  r.  can  rewrite  Equation  IS  as: 

(o  (r) !  ‘  =  R  (  r.r  )  -  <  I  (r)  >‘  O0) 

ri  n  ’ 

It  can  be  shown  (Goodman.  19S5)  that  another  way  to  write  Equation  19.  with 
r ]  =  f;  =  r.  is  given  by: 

Rn  (  r.r  )  =  <  hr)  >  <  Hr)  >  +  |  R  (  r,r  )  |‘ 

where  R  (r.r)  is  the  ATC  of  the  image  amplitude  and.  from  the  relationship  be¬ 
tween  Equation  12  and  Equation  19.  is  equal  to  <  g(r)  (g(r))  *  >  (since  amplitude 
can  be  complex,  in  eeneral).  For  t(r)  and  h(r).  though,  this  will  always  be  a  real 
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function,  so  I  can  write  Rge  (r,r)  =<  |  g  (r)  |2>=  <l(r)  >.  Equation  21  can  be 
written  as: 

Rn  (  r,r  )  =  <  I(r)  >2  +  <  l(r)  >2  =  2<  I(r)  >2  (22) 

Using  this  result  in  Equation  20  gives  (  on  (r))2  =  <  I (r)  >2,  which,  upon  taking 
the  square  root  of  both  sides  and  dividing  by  the  ensemble  averaged  intensity, 
gives  us  the  same  result  as  Equation  9  for  the  contrast.  So  the  rms  of  the  intensity 
equals  the  mean  intensity  even  for  non-stationarv  objects.  The  difference  here  is 
that,  in  general,  the  variance  is  not  a  constant  but  is  a  result  of  the  ensemble 
average  and  the  image  coordinates.  Consider  the  mean  intensity  of  the  image  given 
as.  by  definition.  <  l(r)  >  =  <  g(r)  g(r)  *  >.  However,  as  shown  above,  this  is 
nothing  more  than  the  ATC  of  g(r),  given  as  <  I(r)  >  =  Rg!:  (r.r)  . 

Another  important  property  is  the  ATC  of  the  image.  Returning  to  my  original 
notation,  where  dir)  is  assumed  stationary,  its  ATC  can  be  written  as 

R,  .(  r,.  r„  )  =  R.  .  iT.  -iV)  =  &  (r.-r.)  since  I  assume  the  ATC  peak  to  be  verv 

narrow  for  the  wideband  noise  (6  represents  the  Dirac  delta  function).  The  ATC  of 
the  speckled  object  fin  then  becomes: 

R<;  ( r. .  r:  )  =  <  t  (r  j )  d  (r  j )  t  *  (r  2  )  d  *  (r  :  )  > 

R, i  ( r j .  r;  )  =  <  1 1 r  j )  t  (r  2  )  *  >  R  dd  O' i  -  )  (-?> 

R. ;  < r j .  r;  )  =  t  (r ,  1  t  (r;  )*  8  (r  j  -  r:  ) 
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Note  that  R{f  (rp r2)  is  still  a  function  of  r1  and  r2  and  not  just  ri~r2>  even 
though  d(r)  is  stationary,  since  1  have  retained  the  nonstationarity  of  t(r).  So  the 
image  amplitude  ATC  can  be  expressed  in  integral  form  as: 

Rgc  (ri*r2  >  =  /kro)|2h(ri  -  ro)h(r2  -ro)*dro  (24) 

n0 

where  the  integration  is  over  the  object  plane  IIq  and  a  linear  filter  operation  is 
assumed  with  the  impulse  response  h(r).  Again  letting  rj  =  r2  =r,  and  using  the 
relation  from  a  previous  paragraph  of  <  I  (r)  >  =  Rcc(r,r)  in  the  expression 
above  (Equation  24)  I  can  write  (*  denotes  convolution): 

<  1  (r)  >  =  |t(r)|:  *  |  h(r)  |"  (25) 

From  Goodman  (196$)  this  equation  is  of  the  form  of  the  transfer  function  for  a 
system  illuminated  by  incoherent  light  (where  |t(r)|:  is  the  source  intensity  and 
h(r)  is  the  system  impulse  response). 

Thus,  the  transfer  function  that  determines  the  mean  intensity  in  the  image  is  the 
incoherent  transfer  function  of  the  system,  even  for  an  optical  system  illuminated 
by  a  coherent  diffuse  object  field!  One  can  see  from  Equation  25  that  the  PSD  of  a 
speckle  pattern  has  the  shape  of  the  ATC  of  the  squared  modulus  of  the  pupil 
function  (Goodman.l9$5). 

Labeyrie  and  Dainty  (1973)  use  the  principle  of  Equation  25  to  analyze  atmos¬ 
pheric  speckle.  1  mention  this  because  it  relates  a  practical  application  of  this  ex¬ 
pression.  By  measuring  the  image  intensities  and  summing  their  individual  Fourier 
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transforms  from  each  realization,  and  then  dividing  this  result  by  the  time-aver¬ 
aged  modulus  of  each  transfer  function,  the  Fourier  magnitude  of  the  object’s 
intensity  can  be  reconstructed. 

Miller,  et  al.  (1975)  use  a  more  rigorous  approach  to  come  to  the  same  conclusions 
as  given  above.  They  find  that  the  intensity  is,  in  general,  a  nonstationary  random 
process  so  the  ATC  must  be  a  function  of  two  variables.  They  make  the  point  that 
these  second-order  results  are  good  for  both  Fresnel  (so-called  near  field)  and 
Fraunhofer  diffraction  but  that  the  usual  ATC/PSD  Fourier  relationship  does  not 
necessarih  hold.  For  nonstationary  images,  the  PSD  of  the  intensity  is  given  as: 

OO 

S|(t)-  <l/l  <r)  exp  [-j2irjr]  dr  |-  >  (26) 

_  OO 

where  Sj  (  C ),  the  PSD.  is  the  average  value  of  the  modulus  squared  Fourier-trans¬ 
formed  image  intensity  and  i  is  a  spatial  frequency  vector. 

This  relationship  was  extended  by  April  and  Lowenthal  (1984).  They  showed  that 
Equation  26  is  really  only  one  of  two  ports  of  the  PSD  of  the  image  intensity.  The 
complete  PSD  is  given  as  Sj(C)*=  Sj(£)j  +  S  j  ( t )  2  (ecl-  26  is  s  j  (  fe  )j  )  .  This  is 
derived  simply  by  taking  the  PSD  of  Equation  21.  The  first  term  is  the  object  and 
the  second  term  is  the  speckle  noise.  Sj  ( C  )2  is  written  as: 

or 

S|(  c  );  =ff  |<g  (>']  )g  * (r 2  ) > |  exp  [-j2 nr £(  r  j  -r2  )  ]  dr  j  dr2  (27) 
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The  key  point  in  all  of  this  is  that  the  noise  spatial  frequencies  are  across  the  same 
frequency  region  as  the  incoherent  transfer  function  relationship  given  by  Equation 
26.  This  was  proven  by  April  and  Lowenthal  experimentally.  Using  the  speckled 
image  of  a  Ronchi  ruling,  they  showed  that  it  is  th_  incoherent  cutoff  frequency 
(which  is  twice  the  coherent  cutoff  frequency  (Goodman  (1968)  and  Gaskill))  that 
determines  the  distribution  of  the  spatial  frequencies  present  in  the  diffuse  coher¬ 
ent  light.  Note  that  all  indication  of  the  object’s  presence  is  still  lost  in  a  speckled 
image  if  the  incoherent  cutoff  frequency  is  less  than  the  fundamental  frequency  of 
the  Ronchi  ruling.  Kozma  and  Christensen  (1976)  also  experimentally  verified  that 
the  aperture  of  a  coherently  illuminated  system  must  be  twice  as  large  as  an 
equivalent  incoherently  illuminated  system  to  achieve  the  same  resolution  for  each. 
They  compared  bar  targets  of  varying  spatial  frequencies  illuminated  by  both  inco¬ 
herent  and  coherent  sources. 

To  summarize  all  of  these  points,  it  was  shown  that,  even  using  nonstationary 
statistics,  the  rms  intensity  still  equals  the  mean  intensity  and  so.  for  fully  devel¬ 
oped  speckle,  the  contrast  is  one.  Also,  nonstationary  statistics  do  not  change  the 
fact  that  the  rough  surface  is  a  complex  Gaussian  random  variable  leading  to  a 
negative  exponential  distribution  in  intensity.  Finally,  the  signal  was  shown  to  be 
the  object  intensity  attenuated  by  the  incoherent  transfer  function  of  the  system. 
Therefore,  information  transmission  with  coherent  diffuse  illumination  is  depend¬ 
ent  on  the  incoherent  bandwidth  of  the  optical  system. 


3.5  Multiplicative  Noise  Model 


These  results  are  the  basis  for  the  noise  model  used  by  most  researchers  (Kuan,  et 
al,  1987;  April  and  Arsenault,  1976,  1984,  1986;  Guenther,  et  al,  1978;  Jain  and 
Christensen,  1980).  They  describe  speckle  noise  as  multiplicative  in  the  sense  that 
the  speckled  object’s  intensity  is  a  product  of  the  incoherent  signal  and  a  noise 
function. 

Therefore,  the  amount  of  speckle  noise  is  proportional  to  the  image.  The  signal- 
dependent  nature  of  speckle  is  due  to  the  fact  that  the  rms  amplitude  of  the  speckle 
noise  is  proportional  to  the  image.  This  model  is  described  by  (Tur,  et  al,  1982)  as: 

I,(r)  -  <*I0(r)  ]N(r)  (28> 


where: 

Is  O')  is  the  speckled  intensity 

I«(r)  is  the  incoherent  image  intensity  of  t(r) 

I\(r)  is  the  speckle  noise  intensity  of  d(r) 

cv  is  a  proportionality  constant  dependent  upon  the  system  parameters. 

Tite  above  model  breaks  down  if  the  object  has  details  beyond  the  resolution  capa¬ 
bilities  of  the  coherent  optical  system  used  to  image  the  object.  In  other  words,  if 
the  bandwidth  of  the  image  is  greater  than  the  bandwidth  of  the  optical  system,  the 
multiplicative  model  will  not  be  applicable.  Lim  and  Nawab  (1981)  ensure  this 
condition  will  not  occur  by  sampling  the  speckled  image  coarsely  enough  so  that  all 
points  in  the  degraded  image  will  be  considered  independent.  The  sampling  dis¬ 
tance  must  be  larger  than  the  correlation  length  of  the  speckles  at  neighboring 
pixels  to  ensure  a  spatially  uncorrelated  random  field  in  the  image  (Sadjadi.  1990). 
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This  is  the  case  for  most  speckle  computer  models.  In  this  thesis,  I  will  assume  that 
my  coherent  system  can  resolve  all  of  the  details  of  the  object  surface,  t(r). 

3.6  Speckle  Computer  Models 

Having  resolved  that  a  speckled  image  is  a  signal-dependent,  multiplicative  noise 
process.  1  now  turn  to  how  I  want  to  model  this  on  the  computer.  There  are  several 
methods  I  can  choose  to  do  this.  Before  I  describe  them,  I  will  review  speckle. 

A  speckled  object  results  when  a  coherent,  electromagnetic  wave  reflects  from  a 
rough  surface  which  contains  a  very  large  number  of  scatterers.  This  rough  surface 
varies  randomly  and  uniformly  at  depths  on  the  order  of  a  wavelength  of  the  inci¬ 
dent  radiation.  The  reflected,  dephased,  yet  still  coherent  wave  interferes  with  it¬ 
self  and  speckle  is  detected.  The  object  structure  modulates  the  speckle  and  sur¬ 
face  roughness  greater  than  one  wavelength  of  the  incident  radiation  produces  no 
greater  amount  of  speckle.  Once  the  surface  roughness  exceeds  the  coherence 
length  of  the  source,  though,  speckle  will  cease  since  it  no  longer  has  the  modulo 
2tt  (or  integer  multiples  thereof)  phase  shift.  The  speckles  ‘'wash  out"  due  to  the 
loss  of  coherence  of  the  source.  This  is  the  idea  behind  speckle  reduction  tech¬ 
niques  that  are  performed  physically. 

The  negative  exponential  probability  distribution  of  speckle  intensity  has  been  ex¬ 
perimentally  verified  (Dainty.  1976).  This  probability  density  function  (pdf)  can 
result  from  a  circular,  complex  Gaussian  random  variable  with  the  real  and  imagi¬ 
nary  components  having  zero  means  and  equal  variances  (Goodman.  19S5).  The 


standard  deviation  and  the  mean  of  the  intensity  are  equal  for  fully  developed 
speckle,  giving  a  contrast  (C)  ratio  (an  inverse  SNR)  of  one. 

The  most  direct  way  to  model  this  is  given  by  Guenther,  et  al,  (1978).  The  speckle 
field  can  be  described  as  (using  m,n  as  pixel  locations): 

a  (m,n)  =  aR(m,n)  +  jaj(m,n)  (29) 

where  aR(m.n)  and  a  ]  (m.n)  are  zero  mean.  Gaussian  random  variables 
for  the  real  and  imaginary  components,  respectively,  with  variances  of 
°a“  =  E  { “R"  (m.n)}=  E{aj-  (m.n)}  >  The  speckle  intensity  becomes: 

g  (m.n)  =  aR:  (m.n)  +  ai:(m,n)  =  |a  (m,n)  |:  (30) 

To  model  digital  speckle.  1  first  generate,  for  each  pixel,  a  pseudo-random  pair  of 
Gaussian  random  variables  of  zero  mean  and  equal  variances  based  on  the  field 
strength  (square  root  of  the  intensity)  of  the  object  at  each  pixel.  The  Gaussian 
random  variables  of  zero  mean  and  unity  variance  are  created  by  a  routine  from 
Forsythe.  Figure  1  i-  a  typical  histogram  of  the  Gaussian  random  variable  1  have 
created  for  a  given  seed  (7654321).  I  then  take  the  complex  absolute  value  squared 
of  this  quantity  to  obtain  the  speckle  intensity.  Different  realizations  are  provided 
by  using  different  starting  random  seeds  for  each  new  image.  The  optical  system  is 
ignored  in  this  approach  and  completely  uncorrelated  speckle  is  produced.  1  call 
this  technique  the  Direct  Digital  Model  (DDM).  Its  flowchart  is  shown  in  Figure  2. 


Figure  1.  Typical  Histogram  of  a  Gaussian  Random  Number  Function 

Generated  by  Forsythe  Routine 
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DIRECT  DIGITAL  MODEL 


ajm.n)  =  G  AUS_R  AN2(G,)j 


a(m,n)  =  aRe(m,n)  +  j  a,m(rn,n) 


SPECKLE  INTENSITY:  g(m,n)  r 


a(m,n)  p 


Signal  dependency  occurs  because  the  rms  of  the  field  is 
proportional  to  the  input  signal. 

The  function,  a(m,n),  is  a  circular,  complex,  Gaussian 
random  variable  with  zero  mean  real  and  imaginary 
components  having  independent,  identically  distributed 
variances. 


Figure  2.  Flow  Diagram  for  Direct  Digital  Model 


This  result  can  be  extended  to  include  the  optical  system  by  propagating  the  speck¬ 
le  field  to  the  farfield  and  multiplying  this  by  the  coherent  transfer  function  of  the 
aperture.  Inverse  Fourier  transforming  this  function  (focusing  through  a  lens)  and 
then  taking  the  magnitude  squared  (square  law  detector)  gives  a  speckled  object 
that  includes  the  coherent  point  spread  function  (psf)  of  the  optical  system  (Enloe, 
1967;  Goodman,  1968).  I  call  this  technique  the  Transfer  Function  Model  (TFM). 
Its  flowchart  is  seen  in  Figure  3. 

A  technique  used  by  Marathay,  et  al,  (1989)  and  Voelz,  et  al,  (1988)  is  to  create 
a  random,  uniform  phase  of  modulo  2tt  and  then,  using  Euler's  relationship 
exp[j0]  =  cos0  +  jsin0,  multiply  the  phase  by  the  object  field  at  each  pixel  location. 
From  Friedan,  this  phase  is  given  by: 

f(m.n)  =  tt  ( 2  R — 1 )  (31) 

where  0  <  R  <  1  is  a  uniformly  distributed,  pseudo-random  number  produced  by  a 
random  number  generator.  The  key  here  is  that  the  phase  is  statistically  independ¬ 
ent  of  the  original  object.  1  call  this  technique  the  Random  Phase  Model  (RPM). 
See  Figure  4  for  its  flowchart. 
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TRANSFER  FUNCTION  MODEL 


!C''1(B(i,j)}  =  b(m,n)j 

_ U _ 

SPECKLE  INTENSITY:  g(m,n)  =  I  b(m,n)  f 


H(i,j)  is  the  coherent  transfer  function  of  the  optical 
system 


Figure  3.  Flow  Diagram  for  Transfer  Function  Model 


RANDOM  PHASE  MODEL 


The  phase  is  statistically  independent  of  the  original 
object  and  is  modulo  27t  ( -7t  <  <(>  <  n  ) 


Figure  4.  Flow  Diagram  for  Random  Phase  Model 
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Another  way  to  obtain  g(m,n)  above  is  given  by  Kuan,  et  al,  (1987).  They  describe 
the  complex  quantity  b(m,n)  as: 

b  ( m,  n)  =  2  X  h  (m-i ,  n  -  k  )yj  f(i,k)  exp  [  j  0  (i ,  k) )  (32) 

i  k 

where  h(i,k)  is  the  coherent  psf  of  the  system  for  pixels  i,k.  Kuan  terms  Equation 
32  as  the  single  phase  speckle  model.  He  extends  this  model  to  what  he  calls  the 
multiple  phase  model  by  rewriting  Equation  32  as: 

b(m,n)  =  2Sh(m-i,  n-k)^l(i,k)  a(i.k) 

i  k 

where  a(i.k)  is  a  circular,  complex,  Gaussian  random  variable.  Kuan  prefers  this 
latter  model  because  it  reduces  the  sampling  rate  to  that  which  accurately  reflect? 
the  speckle  statistics  and  ‘"abstracts  the  dependence  of  the  speckle  model  on  the 
object  surfaced'  1  have  modeled  this  by  creating  a(i,k)  as  discussed  for  the  DDM, 
multiplying  this  by  the  object  field  and  then  convolving  this  quantity  with  h(i.k)  in 
the  spatial  domain.  Thus,  h(i.k)  is  a  3X3  spatial  window  of  equal  strengths  (1/9) 
that  moves  across  the  image  correlating  the  object  from  pixel-to-pixel  (see 
Gonzalez  and  Wintzi.  Adjacent  pixels  are  no  longer  uncorrelated.  Notice  that  the 
TFM  uses  a  coherent  transfer  function  that  operates  in  the  frequency  domain  upon 
the  entire  image  simultaneously.  Its  correlation  is  much  weaker  than  this  tech¬ 
nique.  The  magnitude  squared  of  this  quantity  gives  the  speckled  object  intensity.  I 
call  this  model  the  Correlated  Window  Model  (CWM).  Figure  5  shows  the  CWM 
flowchart. 
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CORRELATED  WINDOW  MODEL 


|a(m,n)  =  aRc(m,n)  +  j  ^(m.n); 

U _ 

b(m,n)  =  nh(m-i,n-k)a(i,k)! 

I  k  { 

_ i _ 

SPECKLE  INTENSITY:  g(m,n)  =  i  b(m,n)  J\ 


The  function,  h(i,k)  is  a  3X3  spatial  window  of  equal 
strengths  (1/9). 


Figure  5.  Flow  Diagram  for  Correlated  Window  Model 
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The  CW'M  is  probably  the  most  physically  accurate  of  the  speckle  models  since  the 
assumption  of  the  existence  of  totally  independent  scatterers  on  the  diffuser,  d(r), 
is  not  realistic  for  high  quality  optical  systems  (Fujii,  et  al,  1976).  These  systems 
tend  to  correlate  the  scattering  spots  because  they  are  so  smooth  across  the  sur¬ 
face.  Most  researchers  ignore  this  model  because  it  severely  complicates  the  statis¬ 
tical  nature  of  the  resultant  image  thus  complicating  any  statistical  filters  used  to 
reduce  the  speckle  noise.  As  Kuan  shows  in  his  paper,  the  reconstructed  image  is 
more  blurred  than  the  original,  and  1  found  that  to  be  the  case  in  this  paper,  also. 

As  proven  in  earlier  sections,  a  speckled  object  has  the  following  relationship  (in 
Kuan’s  notation): 

g  (m,n)  =  I(m,n)u(m,n)  (34) 

where  u(m.n)  is  a  signal  independent,  white  noise  process  with  a  negative  expo¬ 
nential  pdf  and  I(m.n)  is  the  incoherent  image  of  the  original  image.  This  leads  to 
another  method  of  generating  speckle  images  that  are  “correct"  statistically.  First 
create  u(m.n),  a  function  that  has  an  intensity  that  obeys  a  negative  exponential 
pdf  with  unity  mean  and  unity  variance.  Then  multiply  this  by  the  original  object 
thereby  creating  a  signal-dependent,  multiplicative  noise  process  with  the  requisite 
negative  exponential  pdf.  1  call  this  model  the  Negative  Exponential  Model  (NEM). 
The  NEM  flowchart  is  seen  in  Figure  6. 


EGATIVE  EXPONENTIAL  PDF  MODEL 


SPECKLE  INTENSITY:  g(m,n)  =  t(m,n)u(m,nj 


The  function,  u(m,n),  is  found  from  the  probability 
transformation  of  a  uniform  pdf  to  a  negative  exponential 
pdf. 


Figure  6.  Flow  Diagram  for  Negative  Exponential  pdf  Model 


The  value  of  u(m.n)  can  be  found  from  the  following  equation  (see  Friedan): 
u(m.nj  =  -  In  (  1  -  R  )  (35 

where  R  has  the  same  meaning  as  in  Equation  31  and  In  is  the  natural  logarithm 


Similar  to  the  development  of  the  NEM  is  the  following  technique.  Arsenault  and 
April  (1976)  showed  that  the  natural  logarithm  of  a  speckled  object  has  an  ap¬ 
proximately  Gaussian  pdf  and  is  additive.  To  take  advantage  of  this  property,  1 
create  a  Gaussian  random  variable,  a(m,n),  as  above  and  add  it  to  the  natural 
logarithm  of  the  original  object,  I(m,n).  I  then  exponentiate  this  sum  to  obtain  the 
speckle  intensity.  To  see  this: 

p(m,n)  =  ln(  I(m,n))  +  a(m,n)  (36) 

by  exponentiation: 

g  (m.n)  =  exp  [p(m,n)]  (37) 

I  call  this  model  the  Additive  Gaussian  Model (AGM)  and  its  flowchart  is  seen  in 
Figure  7a.  Figure  7b  is  the  histogram  of  the  function,  a(m.n),  used  to  create  this 
speckle.  These  last  two  methods  are  statistically  correct  but  they  may  be  difficult  to 
model  physically. 

All  six  of  these  speckle  models  are  imaged,  i.e.,  the  speckled  object  exists  in  the 
image  plane  and  not  the  pupil,  or  Fourier  plane.  Several  authors  (Marathay.  et  al. 
1989;  Voelz,  et  al,  19SS:  Idell.  et  al.  1987)  model  this  latter  speckle  condition  but  it 
is  not  considered  in  this  thesis.  Welford  (1976)  reviews  the  difference  between 
imaged  and  non-imaged  speckle  statistics  and  states  that  similar  statistics  will  re¬ 
sult  as  long  as  the  phase  excursion  is  large,  the  number  of  effective  scatterers  is 
large,  and  the  imaging  system  is  not  teleeentrie  (i.e.,  the  optical  path  length  from 
object  plane  to  image  plane  is  the  same  along  all  principal  rays).  All  optical  sys¬ 
tems  considered  in  this  thesis  are  aberration-free. 


Forces  the  logarithm  of  the  speckle  Intensity  to  have  a 
Gaussian  noise  pdf. 

Figure  7a.  Flow  Diagram  of  Additive  Gaussian  Model 
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Figure  “b.  Histogram  of  Zero  Mean,  Unitary  Variance  Gaussian  Function 
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Dainty  (1980)  lists  several  other  types  of  speckle  that  are  seen  but  these  cases  have 
different  statistics  than  the  modeled  versions  in  this  thesis  and  are  merely  listed 
here: 


1.  Non-circular  Gaussian  speckle  -  the  phase  is  not  uniformly  distributed  from 
—77  to  77  but  the  number  of  scatterers  (N)  is  still  large. 

2.  Small  -N  speckle  -  the  number  of  scatterers  are  small  and  thus  non-Gaussian 
statistics  result. 

3.  White  light  speckle  -  spatially  coherent  yet  partially  temporally  coherent 
radiation  is  used  for  illumination. 

4.  Spatially  partially  coherent  speckle  -  similar  to  stellar  speckle  where  atmos¬ 
pheric  turbulence  provides  the  ‘'speckled"  appearance. 

5.  Depolarized  speckle  -  the  object  depolarizes  the  incident  radiation  upon  reflec¬ 
tion. 

The  next  section  will  review  the  checks  necessary  to  verify  the  speckle  models  and. 
give  examples  of  each. 


3.7  Speckle  Verification 


There  are  four  key  ways  to  ensure  an  image  is  “correctly”  speckled,  i.e.,  it  meets 
the  proper  first  order  statistics  that  define  speckle  noise.  These  four  checks  for 
fully  developed  speckled  objects  are  that  its: 


1.  Field  modulus  histogram  has  a  Rayleigh  pdf. 

2.  Intensity  histogram  has  a  negative  exponential  pdf. 

3.  Natural  logarithm  of  the  speckle  intensity  has  an  approximately  Gaussian 
pdf.  i.e.,  exp[-g-exp(-g)]  pdf  seen  in  Figure  8. 

4.  Contrast  ratio  is  approximately  one. 

Checks  one  through  three  are  shown  in  the  following  figures.  Check  four  can  be 
found  one  of  two  ways.  If  my  image  fills  the  entire  array  1  can  compute  the  global 
mean  and  global  variance  to  obtain  the  global  contrast.  If  not,  1  have  to  obtain 
these  values  locally  by.  for  the  mean: 

i  i 

fj.fm.nj  =  —  X  E  c(  m  +  i.  n  +  j  )  (38: 

9  i=-i  i=-i 

for  the  standard  deviation: 


c  (m.n ) 


i  ,  . 

E  (g(  m  +  i,  n  +  j  ) )‘  -  p(m.n) 
,i=-l  9 


a 


for  the  speckle  contrast  (Crimmins,  19S5): 

1  Nv'  1  o  (m.n) 
C  “  <  N  >2  t:  u(m.n) 


Ac. 


i 


(4  0 ) 
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I  compute  these  values  for  each  speckled  image.  To  compare  each  of  these  model¬ 
ing  techniques  I  use  the  same  truth  object.  For  pure  speckle  statistics,  I  use  a 
uniform  “wall”  as  an  object  that  is  the  same  size  as  the  array  (64X64)  and  has  an 

uniform  pixel  intensity  value  of  one.  For  a  constrained  object  1  use  a  sine2 

(x/16,y/16)  in  a  64X64  array  (see  Figures  9a  and  9b).  Several  different  plots  will 
be  shown  for  each  method,  depending  on  the  parameters  I  want  to  use  to  describe 
the  statistics  of  that  particular  method.  An  explanation  will  accompany  each  plot. 
The  same  starting  random  number  seed  was  given  for  all  models  (seed=7654321). 

Figures  10a  through  lOe  demonstrate  the  DDM.  Figures  lOa-d  show  the  speckle 
statistics  for  a  uniform  “wall.”  The  contrast  is  0.9942375  and  the  skewness  of  the 
logarithm  of  the  intensity  is  1.2944.  Notice  that  the  histograms  match  the  exact 
statistical  criteria  that  represent  their  given  phenomenon  (intensity,  field,  log(inten- 
sity)).  The  speckle  phase  is  a  white  noise  process.  Figure  lOe  is  an  example  of  a 
sinc2(r/16)  speckled  via  the  DDM.  Its  imaged  contrast  is  0.9873512. 

Figures  11a  through  Ilf  demonstrate  the  CWM.  Figure  11a  is  a  depiction  of  a 
speckle  pattern  seen  due  to  a  uniform  “wall.”  Its  contrast  is  1.0278078.  Notice  that 
there  is  no  difference  between  the  statistical  histograms  or  the  phase  of  this 
method  and  that  of  the  DDM.  That  is  because  this  technique  is  similar  to  the  DDM 
except  for  the  3X3  spatial  correlating  window.  The  correlating  makes  the  logarithm 
of  the  speckle  intensity  slightly  more  symmetrical  than  the  DDM  (skew- 
ness=1.2122,  so  closer  to  zero  than  the  DDM  for  the  same  implementation)  indicat¬ 
ing  the  localized  Gaussian  psf’s  are  playing  a  small  role. 


HISTOGRAM  Of  SPECKLE  INTENSITY 
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HISTOGRAM  OF  SPECKLE  FIELD  MODULUS 
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CORRELATED  WINDOW  MODEL 
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Figure  11a.  A  Uniform  “Wall"  Speckled  Via  the  CWM 
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histogram  of  speckle  intensity 


HISTOGRAM  OF  SPECKLE  FIELD  MODULUS 


CORRELATED  WINDOW  SPECKLE  METHOD 


Figures  12a  through  12c  demonstrate  the  TFM.  The  TFM  is  another  modification 
to  the  DDM  using  a  coherent  transfer  function.  The  statistics  of  this  technique  are 
too  dependent  on  the  choice  of  the  transfer  function  to  be  universally  applicable. 
For  example,  Figure  12a  is  the  speckled  sine2  (C=0.9634319)  using  a  rectangular 
aperture  of  width  N=64  pixels.  Figure  12b  is  the  same  object  but  with  a  circular 
aperture  of  diameter  N=64  pixels.  Notice  how  the  speckle  is  much  smoother.  Fig¬ 
ure  12c  is  the  image  if  1  used  an  incoherent  transfer  function  and  the  circular 
aperture  with  the  TFM. 

Figures  1 3a  through  13g  demonstrate  the  RPM.  Figure  13a  shows  the  imaged, 
speckled  sine2  (C=0. 7850590).  Using  the  modification  in  Figure  13b,  where  I 
propagate  the  speckle  field  to  the  farfield  (via  Fourier  transformation  (Goodman, 
1968)).  I  obtain  the  speckle  pattern  shown  in  Figure  13c  (C=1.001926).  Comparing 
these  figures  with  their  contrast  ratios,  it  is  obvious  that  the  best  way  to  use  this 
technique  is  for  pupil-plane  speckle  and  not  for  image-plane  speckle  (see  Yoelz. 
et  a!.  1988  and  Marathay.  et  al.  19S9).  The  skewness  for  the  farfield  case  was 
1.20904.  This  technique  could  also  be  used  to  create  the  speckle  field  and  then  be 
imaged  via  a  coherent  transfer  function  similar  to  the  TFM. 

Figures  14a  through  14e  demonstrate  the  NENl.  The  imaged  speckle  contrast  for 
this  model  is  1.071339.  By  forcing  the  object  to  have  a  negative  exponential  pdf 
with  unit  mean  and  variance,  it  is  not  surprising  that  the  histograms  fulfill  the 
statistical  criteria  so  well.  Figure  14c  is  the  histogram  of  the  speckle  intensity  for 
the  "wall"  object  (C=. 9949181 ).  The  skewness  is  1.3722. 
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RANDOM  PHASE  MODEL 
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g(m,n)  =  |  A(m,n)  I2 
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Fiow  Diacram  for  the  Random  Phase  Model 
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Figure  1  ?c.  Truth  Object  Speckled  Via  the  RPM  in  the  Pupil  Plane 
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HISTOGRAM  OF  SPECKLE  INTENSITY 


HISTOGRAM  OF  SPECKLE  FIELD  MODULUS 
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Figure  13d-g.  Statistical  Histograms  and  Phase  for  the  RPM 
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Figure  14b-e.  Statistical  Histograms  and  Phase  for  the  MEM 


Figures  15a  through  15h  demonstrate  the  AGM.  Figure  15a  is  the  AGM  using  a 
“wall”  for  an  object  (C=l. 214454)  and  Figures  15b-e  exhibit  its  statistics.  Notice 
that  the  logarithm  of  the  speckle  intensity  has  a  Gaussian  histogram  (skew- 

ness=7.562X10'4  as  it  should  given  that  the  object  had  a  completely  singular  inten¬ 
sity  (ln(l)=0).  Figure  15f  shows  the  speckled  sine2  using  the  AGM  (C=l. 102836). 
Figure  15g  show  its  logarithm  of  the  speckle  intensity  histogram  and  Figure  15h 
shows  the  same  statistics  but  for  the  logarithm  of  the  object  intensity  subtracted 

first  (skewness=6.56XlO~4 ).  This  simply  verifies  that  my  model  was  implemented 
correctly  (see  Figure  7a). 

From  these  results,  it  is  clear  that  a  coherently  speckled  object  can  be  modeled  in 
either  the  pupil  or  image  plane.  For  my  bispectral  analysis,  I  will  use  the  DDM. 
CWM,  and  the  AGM  for  the  speckle  modeling.  As  will  be  shown  in  the  2-D  recon¬ 
struction.  the  bispectrum  works  particularly  well  for  the  AGM  case.  Appendix  A 
has  a  tabular  synopsis  of  all  of  these  speckling  techniques. 


ADDITIVE  GAUSSIAN  MODEL 


Figure  15a.  Speckled  “Wall"  Using  the  AGM 


Histocxux  or  sptcKir  i-teN5ii»  mivtoqoam  or  anmu  riEitt  Kwmus 


A/IT 1  M  V?  U'l«t1U>  trr«  '  Mfftfvt  tAlBtllM  **T* 


( ll ) 


(  C  ) 
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Figure  15g-h.  Truth  Object  Speckled  Via  the  AGM 
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4.0  BISPECTRUM  ESTIMATION 
4.1  Overview 

The  definition  of  the  bispectrum,  as  stated  previously,  is  given  as  the  Fourier  trans¬ 
form  of  the  triple  correlation.  A  more  rigorous  definition  than  this  is  given  by 
Nikias  and  Raghuveer  (1987),  where  the  bispectrum  is  the  Fourier  transform  of  the 
third  order  cumulant  sequence  of  a  random  process.  This  definition  satisfies  a 
broader  range  of  statistical  criteria,  using  cumulants,  than  the  correlation  does.  But 
for  purposes  of  this  thesis.  tha  correlation  definition  will  suffice.  Having  stated  this 
in  words  1  now  relate  its  mathematical  description,  using  the  notation  of  Lohmann 
and  Wirnitzer  (1984).  The  triple  correlation  is  given  as: 

OO 

i(^  (Xj,  x,)  =J i  (x)  i  (x  +  x 2 ) i  (x  +  x2)  dx  (41) 

-CO 

where  i(x)  is  a  1-D  object  (or  a  two  dimensional  (2-D)  object  with  x  being  a  2-D 
vector). 

The  bispectrum  is  then  formally  given  as: 

CO  CO 

I(  '•'(  fj .  f  2  )  -J  Ji  ( x  j -  x2 )  exp  [— 27T j (f j  x  j  +  f2  x2  )  ]  dx]  dx2  (42) 

-  CO  — OO 

This  can  be  written  as  (see  Appendix  B): 

1  5  ( f  j  •  f  2  )  =llfj)  1  (f;)I*(f1  +  f2)  (4-'’ 

The  bispectrum,  then,  becomes  the  triple  product  of  the  object  Fourier  transform 
over  spatial  frequencies  fi  and  f2.  From  Equation  4?  it  can  be  seen  that  the 
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bispectrum  is  generally  complex  and  has  redundant  regions  providing  hexagonal 
support  (Nikias  and  Raghuveer,  1987).  This  is  shown  in  Figure  16  (Bartlelt,  et  al, 
1984).  Note  that,  from  Equation  43,  a  1-D  Fourier  object  results  in  a  2-D 
bispectrum  and  thus,  by  inference,  a  2-D  Fourier  object  will  give  a  four  dimen¬ 
sional  (4-D)  bispectrum.  This  can  be  a  concern  for  large  dimensions  due  to  the 
memory  space  required  in  a  computer.  Therefore,  any  algorithm  which  applies  the 
bispectrum  must  attempt  to  resolve  this  by  sampling  only  portions  of  the 
bispectrum,  since  it  has  redundant  information  built  into  it. 


q 


(i)  I<3)(p,q)  =  I(3)(q.p) 

(ii)  I(3)(p,q)  =  I(3>(p,  -  P  -  q) 

(iii)  I(3>(p,q)  =|l(3)(-p,  -  q)  | 


if  i  (x)  is  real 

Figure  16.  Redundant  Regions  of  the  Bispectrum 
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4.2  Properties 


Signal  processing  with  the  bispectrum  has  several  important  advantages.  Unlike 
the  power  spectrum  density  (PSD),  Fourier  phase  information  is  retained  in  the 
bispectrum  so  that  unique  image  reconstruction  can  occur.  Fourier  phase  is  lost  in 
the  PSD  due  to  the  complete  spatial  symmetry  of  the  autocorrelation  while  the 
bispectrum’s  support  retains  this  object  phase.  This  is  quite  easy  to  see  from  the 
definition  of  the  phase  of  the  bispectrum  (see  Appendix  B).  This  property  is  also 
useful  in  identifying  whether  a  system  is  maximum,  minimum,  or  mixed  phase 
(Nikias  and  Raghuveer,  1987). 

Another  advantage  of  the  bispectrum  is  that  it  is  shift-invariant,  i.e.,  it  is  insensi¬ 
tive  to  linear,  translational  motion.  This  property  allows  one  to  use  ensemble  aver¬ 
aging  to  reduce  noise  present  in  the  image  without  concern  for  registering  images. 

When  averaging  M  images  together,  the  noise  decreases  as  the  (Jain  and 
Christensen,  1980).  If  these  images  are  not  perfectly  aligned,  the  resultant  noise- 
reduced  image  will  be  blurred  (Gonzalez  and  Wintz).  Since  the  bispectrum  is  in¬ 
sensitive  to  this  misregistration,  ensemble  averaging  can  be  performed  without  this 
image  blur  if  the  averaging  is  performed  in  the  bispectral  domain.  In  this  thesis, 
this  property  will  be  tested  by  randomly  shifting  the  image  from  frame-to-frame 
while  bispectral  averaging  is  being  performed. 
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Another  property  of  the  bispectrum  has  to  do  with  its  third  order  correlation.  If  the 
noise  present  in  the  image  has  a  symmetric  pdf,  then  its  third  order  moment  is 
zero  in  the  average.  This  is  a  property  of  all  odd-order  correlations.  If  the  noise, 
n(x),  is  signal-independent,  additive  and  stationary,  then  the  noisy  image,  j(x),  can 
be  written  as(Lohmann  and  Wirnitzer,  1984): 

j(x)  =  i(x)  +  n(x)  (44) 

Taking  the  ensemble  average  of  the  triple  correlation  of  Equation  44  gives  (Loh- 
mann  and  Wirnitzer,  1984;  Bartelt,  et  al,  1984;  Freeman,  et  al,  1988): 

<j(3)(  x ! ,  x 2 ) >  =  i(3)(  xlt  x2)  +<  n(3) (  xlt  x2)  > 

+  <n(x)>[i(2)(  Xj  )  +  i(2)(  x2)  +  i(2 }(  x2  -  Xj)]  (45) 

+  i  [<n(2)(  X}  )>  +  <n(2)(  x2)>+n(2)(  x2-  x  j )  ] 

where: 

<  >  denotes  the  ensemble  average 

n(i)  (x)  is  the  autocorrelation  (ATC)  of  the  noise 

n  l xi,  x2)  is  the  triple  correlation  of  the  noise 

i  is  the  image  mean 
<n(x)>  is  the  noise  mean 

i(3)  (xi,  x2)  is  the  triple  correlation  of  the  image. 

Three  undesired  terms  exist  that  need  to  be  dissolved.  If  the  noise  is  assumed  to  be 
zero  mean,  all  values  multiplied  by  that  term  go  to  zero.  Thus  the  ATC  of  the 
image  is  of  no  concern.  If  the  noise  has  a  symmetric  pdf,  its  third  order  moment 

goes  to  zero  and  thus  the  n(j)  (xj,  x2)  term  vanishes.  Finally,  the  i  term  can  be 
zero  if  the  image  mean  is  subtracted  from  the  image  before  processing,  or  if  the 
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image  mean  is  zero  initially.  The  terms  multiplied  by  that  value  are  ignored  if  this 
is  the  case.  All  that  is  left  is: 

<j(3)(  Xj,x2)>  =  i(3)(  XJ.X2)  (46) 

From  this  relationship  the  object  Fourier  magnitude  and  phase  can  be  recovered 
(via  the  bispectrum)  and  then  the  original  object  can  be  reconstructed.  Sun- 
daramoorthy,  et  al,  (1990),  proved  this  same  result  in  the  Fourier  domain.  They 
found  that  as  long  as  the  bispectral  axes  were  not  used  in  the  reconstruction, 
unbiased  estimates  would  be  obtained  for  the  original  object  (if  it  was  in  an  addi¬ 
tive  noise  process  that  had  a  zero  bispectrum). 

Two  key  points  should  be  considered  here  before  I  continue.  One  is  that  the  noise 
above  is  additive  and  that  the  noise  reduction  is  performed  in  the  ensemble.  As 
stated  previously,  speckle  noise  is  multiplicative.  The  task  is  to  convert  this  multi¬ 
plicative  process  into  an  additive  one.  To  do  this  a  homomorphic  transformation 
can  be  taken  on  the  speckled  image.  The  transformation  is  called  homomorphic 
because  a  homomorphism  from  the  set  of  multiplicative  numbers  to  an  additive 
group  of  numbers  has  occurred.  A  homomorphic  transformation  is  one  in  which 
the  natural  logarithm  of  the  image  intensity  is  taken  before  any  further  processing 
is  performed.  The  multiplicative  noise  process  then  becomes  additive  due  to  this 
logarithmic  operation.  Several  authors  review  this  technique  in  general  (Oppen- 
heim  and  Schafer,  1975;  Lim,1987;  Arsenault  and  Levesque,  1984;  and  Gonzalez 
and  Wintz.  1987)  while  others  have  reviewed  this  transformation  for  speckle  in 
particular  (Jain  and  Christensen, 1980;  Lim  and  Nawab,  1981;  and  Sadjadi,  1990). 
Arsenault  and  April  (1976)  have  shown  that  an  object  in  speckle  noise  becomes 
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additive  and  signal-independent  with  an  approximately  Gaussian  pdf  when  under¬ 
going  the  transformation  of: 

D  =  -In  (I)  (47) 

where  D  is  the  density  (a  term  left  over  from  the  film  grain  studies)  and  I  is  the 
speckle  intensity.  Therefore,  theoretically  at  least,  performing  a  homomorphic 
transformation  on  Equation  34  will  give  me  the  result  of  Equation  44.  After  the 
homomorphic  transformation  is  performed  for  each  speckled  image,  these  images 
should  be  ensemble  averaged.  Guenther,  et  al,  (1978)  shows  that  this  should  re¬ 
duce  the  noise  level  by  the  yJTi,  where  N  is  the  number  of  images  processed.  How 
this  will  be  done  using  the  bispectrum  is  discussed  next. 

4.3  Algorithms 

In  this  section  I  will  review  the  algorithms  used  to  reconstruct  multiple,  speckle-de¬ 
graded,  jittered  images.  To  begin  I  refer  to  the  bispectrum  definition  of  Equation 
43,  written  here  as: 

I!3,(f,,f:)  =  |l(3' (fi,f2)|  exp  [Mfi,f2)l  (4S) 

This  function  has  both  magnitude  and  phase,  broken  into  their  respective  compo¬ 
nents  here  as: 

|l(3,(f,.  f; )  I  -  |l(fi)|  |l(f2)l  I  I*  (fi  +  f2)l  («) 

<|((fi,f2)  ■=  4>  (f i )  +  <l>(f2)  -  <t>(fi  +  f2  )  (50) 
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For  a  noisy  sequence  j(k)=i(k)+n(k),  from  0  <  k  <  N-l,  its  Fourier  transform  is 
given  as: 

J(M«  E  j(n)  exp[-j2TTn\]  ^ 

n  =  -  oo 

For  M  images,  the  ensemble-averaged  bispectrum  of  Equation  51  becomes: 

A  1  m 

J<  (52) 

A 

where  j(3)(Xr  X-,)  is  the  bispectrum  estimate  of  the  noisy  image.  The  ultimate 

A 

goal  is  to  obtain  a  bispectrum  estimate  J(3)  that  is  as  close  as  possible  to  the 
image  bispectrum  1^  (  Xlt  X2)  ,  the  Fourier  transform  of  Equation  46.  Once  this 
is  attained  image  recovery  can  be  performed  by  the  following  algorithms. 

4.3.1  One-Dimensional  Algorithm 

For  1-D  images  I  use  a  different  algorithm  each  for  the  magnitude  and  the  phase 
reconstruction.  The  magnitude  algorithm  is  from  Sundaramoorthy,  et  al,  (1990) 
while  the  phase  algorithm  is  from  Bartelt,  et  al,  (1984).  Matsuoka  and  Ulrych 
(1984)  discuss  several  other  phase  algorithms  that  could  be  used  for  1-D  images 
(seismic  wavelets)  but  the  recursive  algorithm  of  Bartelt  appears  to  be  sufficient.  1 
will  first  develop  the  1-D  algorithm  for  a  1-D  sequence  from  0  <  m.n  <  N-l. 
letting  B(m,n)  denote  I(3!(Xlt  X2)  and  H  denote  I.  With  these  substitutions  1  get: 

B(m,n)  =  H(m)H(n)H*(m+n)  (53) 
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where  B(m,n)  is  the  bispectrum  of  H  ( the  Fourier  transform  of  the  1-D  sequence). 
Letting  m=l,  Equation  53  becomes: 

B(l,n)  =  H(l)H(n)H*(l+n)  (54) 

Since  the  object  is  assumed  to  be  real  its  Fourier  transform  is  Hermitian,  i.e., 
f(x)  =  f  *  (— x)  (Gaskill).  Hermitian  means  that  the  complex  function  has  a  real 
component  which  is  even  (e(x)  =  e(-x))  and  an  imaginary  component  which  is  odd 
(o(x)=-o(-x)).  This  implies  I  only  have  to  sample  the  bispectrum  from  Equation  54 

for  1  <  n  <  — — 1  and  then  I  can  later  extend  the  reconstructed  object  from  these 

samples  for  —  to  N-l.  The  value  of  H(0)  is  the  image  mean  (Gonzalez  and 

Wintz).  If  this  value  is  zero,  or  the  mean  is  subtracted  from  the  image  before 
Fourier  transformation,  this  term  could  cause  magnitude  reconstruction  problems, 
since  many  magnitude  algorithms  have  this  term  as  a  divisor  (see  Bartelt,  1984). 


Once  the  bispectrum  values  of  Equation  54  have  been  calculated  and  summed  for 
M  images,  I  can  reconstruct  the  Fourier  magnitude  and  phase  of  the  original  ob¬ 
ject.  For  the  magnitude,  I  simply  write  Equation  54  in  terms  of  its  magnitude, 
letting  n=j-l,  as  (for  real  objects): 


_  B(  1,  j  -1  ) 

H  ( 1)|  |  H(j-l) 


(55) 


where  the  value  of  j  goes  from  2  to  -j-.  For  values  of  j  >  |h(  j)|  =  |  H ( N  —  j ) | . 

Notice  that  I  have  the  I H(  1 ) I  term  as  a  constant  to  determine  before  I  can  recur- 
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sively  find  the  other  values  of  |H(n)|.  This  value  is  given  from  Sundaramoorthy,  et 
al,  as: 


6  |b(i,i)|3  |b (i, 3 ) | 
|H(1)I%  lB(1.2)|  |B(2.2)| 


(56) 


Note  that  an  additional  bispectrum  term  has  arisen  in  this  formulation.  Therefore 
the  constant  B(2,2)  must  be  calculated  and  summed  during  the  calculation  of  the 
other  bispecturm  terms.  From  Equation  53  this  becomes: 


B(  2,  2  )  =  H(2)2H  *  (  4) 


(57) 


From  this  I  can  find  its  magnitude  and  use  its  value  in  the  recursive  routine  (B(2,2) 
is  not  needed  in  the  phase  recursion  routine).  Since  this  technique  requires  only 

samples  of  the  bispectrum  along  the  m=l  line  and  the  B(2,2)  sample  (for  an 

image  of  size  N),  considerable  savings  can  be  had  both  in  memory  allocation  and 
computation  time. 


The  phase  portion  of  the  bispectrum  reconstruction  is  also  recursive.  The  calcu¬ 
lated  values  of  the  bispectrum  phase  are  an  extension  of  Equation  50  used  in 
Equation  54,  written  as: 

\j/(  l,n)  =  <>(1)  +  4>(n)  -  4>  (  1  +  n  )  (5S) 

where,  as  before.  \j/  is  the  bispectrum  phase  and  4>  is  the  object  phase.  Letting 
i  =  n-1,  these  values  are  calculated  only  from  1  <  i  <  I  can  rewrite  Equation  58 
to  determine  the  values  of  4>  as: 
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<|>(i)  =  -^(l,i  -1)  +  4>  (1) 

+  4>  ( i  -1 ) 

(59) 

-M  -i 

where:  4>(1)  =  EM  1,  n) 

n=l 

(60a) 

4>  (i)  =  -4>(N-i )  for  i  > 

?♦> 

(60b) 

Values  for  both  the  magnitude  and  the  phase  for  n  >  -^-+1  to  N-l  can  be  found  by 
taking  advantage  of  the  Hermitian  property  of  the  Fourier  object.  The  value  of 
<K— )  is  set  initially  at  zero  since  it  does  not  affect  the  reconstruction  but  merely 
the  placement  of  the  reconstructed  object  (Matsuoka  and  Ulrych,  1984;  Bartelt, 
1984).  The  final  values  of  J H(  n)  |  and  4>(n)  can  then  be  combined  using  Euler’s 
relationship  and  the  inverse  Fourier  transform  taken  to  obtain  the  reconstructed 
object. 

4.3.2  Two-Dimensional  Algorithm 

The  1-D  algorithm  developed  previously  worked  extremely  well  on  all  types  of 
noise,  as  well  as  speckle  (see  Results  section)  but,  being  a  1-D  algorithm,  it  did 
not  extend  well  to  2-D  images.  By  extension  I  mean  continuing  to  use  this  1-D 
approach  on  2-D  images  by  first  “raster-scanning”,  or  peeling,  the  image  before 
obtaining  its  bispectrum.  If  the  2-D  image  was  noise-free,  this  technique  worked 
perfectly  for  jittered,  image  reconstruction.  All  2-D  objects  of  any  extent,  con¬ 
straint,  or  phase  type  could  be  reconstructed  with  this  method.  Once  a  noise  was 
added,  though,  even  in  small  amounts,  reconstruction  would  not  occur.  Since  noise 
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essentially  creates  phase  errors  in  the  object  phase  (Marron,  et  al,  1990),  these 
phase  errors  are  propagated  and  therefore  increased  when  scanning  an  image  re¬ 
cursively  over  its  NXN  extent.  This  is  because  you  have  given  the  phase  errors  a 
greater  amount  of  data  to  accumulate  in  as  you  scan  the  array  as  opposed  to  a 
1XN  image. 

The  2-D  magnitude  reconstruction  algorithm  I  used  is  from  Dianat  and  Raghuveer 
(1990).  Their  approach  is  developed  as  follows.  Write  Equation  53  for  2-D  images 
as: 

B(i,j;m,n)  =  H(i,j)H(m,n)H*(i+m,j+n)  (61) 

where  the  B(i,j;m,n)  is  the  4-D  bispectrum  of  the  2-D  Fourier  object  H  for  pixels 
m,n,i,j.  To  simplify  to  a  2-D  problem,  let  i  =  m,  j  =  n,  to  get,  only  for  2-D  sub¬ 
plane: 

B(m,n)  =H(m,n)2  H*(  2m,  2n  )  (62) 

Taking  the  natural  logarithm  of  the  magnitude  of  both  sides  of  Equation  61  gives: 

In  [|B(m.n)|  ]  =  2  In  [  |  H(m,n)|]  +  In  [  |  H  (2m,  2n)|]  (63) 

By  describing  the  above  equation  in  terms  of  its  discrete  Fourier  transform  (DFT, 
see  Equation  51)  and  equating  like  coefficients  on  both  sides,  the  following  rela¬ 
tionships  can  be  seen  as  (where  the  ~  term  denotes  the  DFT  of  the  natural  loga¬ 
rithm  of  the  magnitude): 
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1.  For  m,n  odd: 


2. 


|  rs^ 

H(m,n)  =  ~2  B(m,n) 

For  HI  or  —  odd: 

2  2 


(64) 


H(m,n) 


1  ufm  n 

"  t  H  It  .t 
2  ^  2  2 


(65) 


3. 


4. 


Repeat  Step  2  recursively  for  cases  when 
log2(N)-l. 

Finally,  compute: 


is  odd  for  L  =  2,3, 


H(0,0)  =  -i  [B(0,0)-fi(o,-K)  (66) 

Once  I  have  calculated  all  values  of  H(m,n)  recursively  from  the  bispectrum  using 
the  above  relationships,  I  exponentiate  this  value  to  obtain  the  Fourier  magnitude 

|H(rn  n)  |.  I  have  found  this  algorithm  to  work  practically  flawlessly  for  magnitude 
reconstruction.  As  will  be  discussed  in  detail  later,  it  is  the  2-D  phase  that  will 
determine  the  true  shape  of  the  reconstructed  object.  The  above  routine  also  can 
be  extended  to  the  phase  but  it  is  not  as  robust  in  phase  recovery  as  the  magnitude 
algorithm  is.  Thus,  for  phase  reconstruction,  I  used  a  2-D  recursive  routine  similar 
to  the  1-D  in  development  but  with  some  added  calculations. 
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The  two-dimensional,  recursive,  phase  reconstruction  algorithm  is  a  2-D  extension 
of  the  Lii-Rosenblatt  technique  described  by  Matsuoka  and  Ulrych  (1984).  This 
algorithm  requires  the  H(0,1)  and  H(1,0)  subplanes  as  a  minimum  because  these 
subplanes  represent  the  lowest  spatial  frequencies  in  the  2-D  object.  This  being  the 
case,  I  let  i=0,  j=l  in  Equation  61  to  get: 

B(m,n)  =  H(0,l)H(m,n)H*(m,n+l)  (67) 

Similarly,  for  i=l,  j=0  I  get: 

B(m,n)  =  H(l,0)H(m,n)H*(m+l,n)  (68) 

Finally,  I  need  to  calculate  a  single  bispectrum  constant  to  complete  the  calculated 
values  necessary  to  reconstruct  the  image.  This  value  represents  the  limit  of  the 
subplanes  given  above  that  I  will  need  to  sample: 

B  (-2L.0;  0±)  =  H  (-|-,o)  H  (o,-|-  )  H-  (-bL.-SL)  (69) 

To  calculate  the  Fourier  phase  I  use  the  following  relationships,  once  I  have  ob¬ 
tained  the  bispectrum  for  2  <  m,n  <  -j--l: 

1.  Using  \J /  calculated  from  Equation  67  for  m=0,  2  <  n  <  -~1, 

4>(0.n)  =  4>(0.1)  +  4>(0,n-l)  -  \J/(0,n— 1)  (70) 

-i 

where:  4>(1.0)  =  -^  E^(0,n) 

^  m=l 
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4>(0,n+~)  =  -  <K0A-n)  (72) 

2  2 

4.  Using  Equation  71  for  1  <  m,n  <  -—--l, 

4>(m+NV2)  =  -  4>(N72-m,0)  (73) 

5.  Using  y\i  calculated  from  Equation  67  for  1  <  m  <  N-l,  1  <  n  < 

4>(m,n)  =  4^(0. 1)  +  4>(m,n-l)  -  t|/(m>n)  (74a) 

4>(m,n+— )  =  -  4>(N-m.— -n)  (74b) 

2  2 


6.  Finally,  the  following  constant  values,  which  affect  the  position  of  the  object, 
are  set  at: 


4>(0,0)  =  0 
<t>(y,0)  =  0 

4(0,— )  =  0  (or  it  if  B(— ,0-,0.— )  is  less  than  zero). 
2  2  2 
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As  one  can  see  from  this,  finding  the  phase  recursively  for  2-D  images  is  a  bit 
more  involved  than  for  1-D.  Once  both  the  Fourier  magnitude  and  phase  have 
been  determined,  I  combine  them  and  inverse  Fourier  transform  them  to  get  the 
object.  Examples  of  both  1-D  and  2-D  image  reconstruction  will  be  found  in  the 
next  section. 
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5.0  RESULTS 


In  this  section  I  report  my  findings  from  the  computer  simulations  I  have  per¬ 
formed  on  shift-invariant  image  reconstruction  of  speckle-degraded  images  using 
the  bispectrum.  1  have  separated  the  results  into  two  sections:  1)  non-speckled, 
and  2)  speckled.  This  is  because  of  certain  results  I  obtained  from  the  algorithms 
used  above  that  were  totally  unrelated  to  the  speckle  problem. 


In  Appendix  C.  I  have  given  an  overall  flowchart  of  the  2-D  reconstruction  algo¬ 
rithm  used  for  jittered,  speckled  images.  It  is  directly  applicable  to  both  unspeckled 
images  and  all  of  the  1-D  cases.  The  pieces  of  the  algorithm  that  make  up  the 
entire  program  have  been  described  previously. 


> 
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When  I  began  this  thesis  my  first  question  was:  “Is  the  bispectrum  robust  enough 
to  handle  speckled  images?"  The  bispectrum  has  recently  been  used  for  tonal  noise 
and  additive  Gaussian  noise  (Sundaramoorthy,  et  al,  1990)  as  well  as  photon- 
limited  background  noise  with  SNR  <  1  (Newman  and  Van  Yracken,  1989).  The 
bispectrum  has  been  used  in  the  last  few  years  (see  Ayers,  et  al,  198S:  Bartelt.  et 
al.  1984;  Lehmann,  et  al.  1 9S3:  and  Freeman,  et  al,  1988)  to  retrieve  astronomical 
objects  from  atmospheric  speckle  (described  by  Dainty  (1971,  1973)  and  Labeyrie 
(1970)).  Marathay  (1986)  and  Sato,  et  al,  (1978)  have  investigated  methods  to 
physically  perform  correlations  to  help  retrieve  objects.  Thus,  the  recent  interest  in 
bispectrum  estimation  to  reconstruct  objects  in  noise  led  to  this  thesis.  The  concern 
for  coherent  speckle  is  that,  by  definition,  its  phase  is  a  modulo  2tt  random  phase. 


I 


white  noise  process.  Is  this  too  big  of  a  challenge  for  a  phase  retrieval  routine  such 
as  the  bispectrum? 

Another  question  was:  “Can  ensemble  averaging  be  performed  on  multiple, 
shifted,  and  possibly  blurred  images  in  noise  without  image  blurring  in  the  ensem¬ 
ble  occurring  as  a  result?”  This  problem  can  occur  when  imaging  an  object  moving 
against  a  uniform  background  that  is  recorded  by  videotape.  The  same  image  is 
taken  over  many  frames  (at  the  standard  video  rate  of  30  frames  per  second)  and 
trying  to  register  these  perfectly  without  introducing  blurring  is  a  formidable  task. 
Thus,  the  shift-invariant  nature  of  bispectrum  estimation  makes  it  an  excellent 
candidate  to  perform  this  function. 

Therefore,  this  thesis  was  proposed  as  a  conceptual  study  of  image  reconstruction 
of  speckled  images.  Since  ensemble  averaging  of  the  speckled  images  is  the  MLE 
of  the  noise-free  object,  performing  this  in  the  bispectral  domain  while  jittering  the 
object  in  the  spatial  domain  should  test  its  shift-invariant  capabilities  as  well  as  its 
phase-retrieval  capabilities. 

5.1  Non-Spccklctl 

When  I  began  this  task  I  did  not  foresee  that  I  would  have  problems  reconstructing 
a  noise-free  object  and.  in  most  cases.  I  did  not.  There  are  certain  conditions  that 
can  occur  that  can  either  distort  the  reconstructed  object  or  preclude  its  reconstruc¬ 
tion  at  all  and  I  will  discuss  those  conditions  next.  These  problems  were  found  by 
thoroughly  exercising  the  algorithms  for  various  types  of  constrained  objects  in 
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various  extents  and  viewing  the  results.  The  biggest  problem  was  phase  distortion 
of  the  reconstructed  object  and  I  will  review  that  next. 

5.1.1  Phase  Distortion  Problem 

Bispectrum  estimation  results  in  a  linear  phase  shift  being  introduced  into  the  re¬ 
constructed  object,  given  as  H(\)exp[-j27ra\],  where  H(X)  is  the  Fourier  domain 
representation  of  the  reconstructed  object  and  exp[-j277Q;X]  s  the  linear  phase  shift 
(a  is  a  real  constant).  This  phase  shift  can  be  seen  in  the  reconstructed  objects  of 
the  following  figures. 

Figure  17a  is  the  reconstructed  version  of  Figure  9a.  Figure  17b  is  a  2-D  contour 
plot  of  Figure  9b  while  Figure  17c  is  the  recentered  object  using  a  program  1  wrote 
that  centers  the  image  by  wrapping  it  around  given  the  amount  of  shift  desired  in 
both  the  x  and  y  direction.  Notice  that  the  reconstructed  values  go  negative.  This  is 
an  artifact  of  the  2-D  bispectrum  reconstruction.  Figure  18a  is  a  1-D  slice  of  the 
object  in  Figure  26a  and  the  reconstructed  object  is  in  18b.  The  reconstruction 
there  is  exact  except  for  the  phase  shift.  This  phase  shift  acts  as  a  filter  that 
merely  introduces  a  shifting  of  the  object  in  the  spatial  domain.  The  1-D  tribar 
target  of  Figure  19a  did  not  show  a  phase  shift  upon  reconstruction  (Figure  19b). 
probably  because  a  was  zero  so  the  phase  shift  was  zero.  To  explain  this  1  will  use 
a  linear  systems  approach  (Gaskill).  For  an  object.  f(x).  being  convolved  with  a 
filter.  h(x).  I  can  write  (where  (*))  denotes  the  convolution  operator): 

g(x)  =  fix)  *  h(x)  (75) 

and  its  Fourier  transform  relationship  is  (Gaskill): 

G  (  a  )  =  Fi.mHu)  (76) 
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Figure  17a.  Three-Dimensional  Plot  of  Bispectrally-Reconstructed 

Sinc2(X/16,Y/16) 


17b.  Two-Dimensional  Contour 
of  (a) 


17c.  Recentered  Version  of  (b) 


TRUTH  OBJECT:  1-0  SLICE  OF  S/C  (GRYI 


Figure  18a.  One-Dimensional  Truth  Object 


reconstructed  i-a  slice  of  s/c  igryi 


Figure  18b.  Bispectrally  Reconstrucied  Version  of  (a) 
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original  object  (tbibar.  m>S21 
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As  we  have  seen,  each  Fourier  quantity  has  both  magnitude  and  phase  since  the 
Fourier  expression  is,  in  general,  complex.  Thus,  I  can  write  Equation  76  as: 

G(M  =  |f(X)|  |  H(\)|  exp[-j{<J>F(\)  +  4>h(\)}]  (7?) 

But  for  a  phase  filter,  |  H(  \)  |  =  1  and  Equation  76  is  written  as: 

G(\)-F(\)exp[-j4>H(M]  (78) 

If  x)  =  2ttq;X,  then  the  output  function,  g(x)  would  merely  be  a  shifted  ver¬ 
sion  of  the  input,  since,  by  definition,  h(x)  =  8(x  -  a)  in  Equation  75  (where  8  is 
the  Dirac  delta  function  (Papoulis)): 

g(x)  =  f (x)  *  8(x  -  a)  =  f(x  -  a)  (79) 

There  is  no  phase  distortion,  then,  if  the  phase  transfer  function  is  a  linear  func¬ 
tion  of  X.  Also,  if  4>h(\)  is  a  constant  no  phase  distortion  will  occur  (Gaskill). 

For  example,  if  $H(X)  By ,  the  real-valued  input  becomes  an  imaginary-valued 
output.  Similarly,  for  <1>  (\)  =  n2tr  (n  is  an  integer),  g(x)  «  f(x)  and  for 
$„(  X)  =  (  2n  +1)t t,  g(x)  =  -f(x).  The  1-D  tribar  target  of  Figure  18a  did  not  show 

Jri  w  w  w 

a  phase  shift  upon  reconstruction  probably  because  of  the  spatial  frequency  of  the 
tribar  was  such  that  <t>  (X)  =  n2r r. 

The  importance  of  phase  information  cannot  be  overstated  when  reconstructing  an 
object.  In  general,  the  magnitude  gives  you  information  regarding  the  strengths  of 
the  signal  but  the  phase  gives  you  information  about  the  structure  of  the  signal. 
Reconstructing  an  object  merely  from  its  Fourier  magnitude  is  similar  to  taking  the 
autocorrelation  of  the  object  (Goodman,  1985).  This  is  the  result  of  the  square-law 
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detection  aspect  of  image  recording.  As  shown  above,  phase  records  the  location 
of  the  image  whereas  the  magnitude  does  not.  A  simple  example  of  this  is  seen  by 
the  following  three  objects  where  (Dainty  and  Fienup): 

f(x,y)  is  the  original  object 

f(x  -  xo,  y  -  yo)  is  the  object  shifted  by  xq,  xo 

-f(x  -  xo,  y  -  yo  is  the  object  shifted  and  reflected  about  origin. 

All  three  objects  have  the  same  Fourier  modulus  but  all  3  are  not  located  at  the 
same  position  or  orientation.  The  phase  problem  becomes,  then,  as  posed  by  Bates 
and  McDonnell,  does  there  exist  more  than  one  object  belonging  to  the  set  of  all 
possibilities  for  the  phase?  If  the  phase  can  be  retained  exactly,  there  is  no  ambigu¬ 
ity  in  the  object  reconstruction  (Oppenheim,  et  al,  1982). 

An  illustration  that  compares  the  relative  importance  of  phase  versus  magnitude 
for  image  reconstruction  is  given  by  Oppenheim  and  Lim  (1981).  They  provide  the 
example  of  the  reconstruction  of  the  atomic  structure  of  a  L-tyrosine  HCL  crystal 
image  from  its  exact  phase  data.  They  compare  three  succeedingly  worse  cases  for 
the  magnitude:  1)  magnitude  of  one;  2)  magnitude  equal  to  an  average  of  several 
different  L-tyrosine  HCL  atoms;  and  3)  magnitude  from  a  totally  different  atom! 
In  all  cases  the  original  L-tyrosine  HCL  atomic  structure  was  reproduced  from  the 
correct  phase  but  erroneous  magnitude  information. 
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For  further  proof  of  the  importance  of  phase  in  reconstructing  an  object,  I  per¬ 
formed  my  own  experiment  using  the  truth  objects  of  Figure  20a  (“battleship”) 
and  Figure  20b  (“Lena”).  I  combined  the  battleship’s  phase  with  Lena’s  magnitude 
and  obtained  Figure  21a.  Similarly,  I  combined  Lena’s  phase  with  the  battleship’s 
magnitude  and  created  Figure  21b.  Clearly,  the  phase  is  key  when  reconstructing 
an  image. 

Returning  to  the  phase  filter,  if  the  recovered  phase  is  distorted,  the  reconstructed 
object  will  be  distorted.  Phase  distortion  occurs  when  the  “phase-filter”  aspects  of 
the  operation  cause  the  spectral  components  of  the  original  signal  to  be  misplaced 
from  their  original  spacing  (i.e.,  the  amount  of  shift  of  the  spectral  components  is 
not  linear  with  increasing  frequency  (Gaskill)).  When  phase  distortion  occurs,  you 
do  not  have  the  simple,  shifted,  delta  function  convolution  described  in  Equation79 
for  integer  values  of  a.  This  relationship  is  merely  a  special  case  of  the  more 
general  result: 

g(x)  =  f(x)  *  sinc(x-a)  (80) 

where  sinc(x)  =  S'^.^X  ^  (see  Gaskill)  and  a  delta  function  results  at  the  limit  (see 
Papoulis,  p.  280-281).  If  |  ot  |  5.  (Dianat  and  Raghuveer,  1990)  then  nonin¬ 
teger  values  of  the  phase  filter  exist  and  thus  nonlinear  increases  in  the  phase  shift 
occur.  For  values  of  a  >  -£■  up  to  multiple  integer  values,  the  shift  is  merely  a  tt  or 

2tt  phase  shift  and  further  distortion  does  not  occur.  The  reconstructed  object  ap¬ 
pears  to  be  blurred  by  a  convolution  with  a  sine  function.  Examples  of  this  are 
seen  in  Figures  22  through  26. 
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Figure  20a.  “Battleship” 
Used  for  Magnitude/Phase 
Comparison 


Figure  20b.  “Lena” 
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Figure  21a.  Recon¬ 
structed  Object  Using 
“Battleship”  Phase  and 
“Lena”  Magnitude 


Figure  21b.  Recon¬ 
structed  Object  Using 
“Lena”  Phase  and  “Battle¬ 
ship”  Maginitude 
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Figure  22a  is  the  same  tribar  as  in  Figure  19a  except  that  now  it  is  zero-padded. 
Its  reconstructed  image  is  in  Figure  22b.  Notice  that  this  image  exhibits  the  phase 
shifting  and  the  sine  blurring  (which  causes  the  negative-valued  artifacts  seen,  for 
example,  in  Figure  17).  For  the  1-D  case  I  removed  the  sine  blurring  with  data 
windows,  which  I  will  discuss  shortly. 

For  2-D  images,  the  sine  blurring  was  much  worse.  I  tried  several  different  phase 
reconstruction  algorithms  (Dianat  and  Raghuveer,  1990;  Swami  and  Giannakis, 
1988)  before  settling  on  the  one  described  previously  and  it  still  shows  sine  blur¬ 
ring  for  certain  constrained  objects.  For  example,  Figure  23a  is  a  tri(-f-,-^)  truth 
object  reconstructed  in  Figure  23b.  Notice  the  severe  sine  blurring  of  the  object.  I 

now  change  the  constraint  of  the  object  to  be  tri(“,-^)  and  the  reconstructed 

object  of  Figure  24a  is  exact,  save  for  the  introduced  phase  shift.  My  shift-image 
program  then  recenters  the  object  for  viewing  (Figure  24b). 

Not  only  is  object  constraint  a  problem  but  also  the  object  shape  can  cause  sine 
blurring  in  the  reconstructed  image.  Figure  25a  is  what  I  call  a  “binary”  spacecraft 
because  its  gray  levels  change  value  abruptly.  Figure  25b  is  the  reconstructed  ob¬ 
ject.  Notice  that  the  sine  blurring  has  actually  distorted  the  object  and  makes  it 
appear  noisy.  Even  for  smoothly  varying  objects,  such  as  my  “gray”  spacecraft  of 
Figure  26a  (so-calied  because  there  are  several  gray  levels  in  the  image  that  vary 
smoothly),  you  can  see  a  ripple  of  sine  blurring  in  the  reconstructed  object 
(Figure  26b). 
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TRIBAR  IN  ZERO-PADDED  ARRAY  tN-64) 


Figure  22a.  Original  Object 


RECONSTRUCTED  TRIB4R  SHOWING  SINC  BLUR 


Figure  22b.  Bispectrally  Reconstructed  Version  of  (a) 
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Figure  25b.  Bispectrally  Reconstructed  Version  of  (a) 


One  way  to  avoid  the  sine  blurring  in  the  2-D  case  was  to  use  the  1-D  algorithm 
and  raster-scan  the  image  before  calculating  its  bispectrum.  Figures  27a  shows  the 
binary  spacecraft,  its  jittered  representation  is  in  Figure  27b  from  and  the  recon¬ 
structed  object  is  in  Figure  27c.  Notice  that  there  is  no  sine  blurring  of  the  recon¬ 
structed  image  (I  have  recentered  it  after  image  reconstruction).  No  noise  was 
present  so  there  was  no  reconstruction  problems  when  applying  the  1-D  raster- 
scan  method  to  the  2-D  image. 

5.1.2  Windowing  the  Data 

In  the  1-D  case  certain  bispectrum  reconstructions  resulted  not  only  in  sine  blur¬ 
ring  but  also  in  incomplete  reconstructions,  i.e.,  it  would  lead  to  an  indeterminate 
condition  (divide  by  zero).  For  the  1-D  case,  I  believed  these  occurrences  were 
happening  because  of  errors  in  sampling  the  spectrum.  This  was  because  I  noticed 
that  whenever  I  zero-padded  an  array,  the  reconstructed  object  often  had  sine 
blurring  or  it  would  not  reconstruct  at  all.  To  overcome  this,  I  decided  to  use  data 
windows  and  they  worked  as  advertised. 

To  review,  data  windows  are  used  to  relieve  leakage  or  “ringing”  artifacts  when 
spectrum  data  is  aperiodically  sampled.  Data  windows  help  by  forcing  the  sampled 
spectra  to  have  a  period  of  N  extent  (where  N  is  the  length  of  the  data  record)  even 
though  the  object  itself  may  have  an  aperiodic  spectra.  The  formulas  used  for  the 
data  windows  were  taken  from  Stearns  and  David.  In  particular,  I  used  the  Bartlett 
(triangular)  and  the  tapered  rectangular  data  windows. 
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Figure  27a.  Original  Object  Figure  27b.  Representative  Object 

Jittered  Randomly  In  The  Field  Of 
View  For  10  Frames-No  Noise 


a 
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Figure  27c.  Bispectrally  Reconstructed  Version  of  (a) 
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Figure  28  is  an  example  of  using  a  tapered  rectangular  window  on  the  object  in 
Figure  22a  before  calculating  its  bispectrum.  Compare  this  result  to  that  of  Figure 
22b.  This  implies  to  me  that  the  sine  blurring  in  the  1-D  case  is  actually  ringing 
due  to  aperiodic  sampling  and  not  the  same  phenomenon  seen  in  the  2-D  recon¬ 
struction. 

Figure  29a  is  an  example  of  an  object  that  would  not  reconstruct  with  the  1-D 

algorithm  unless  a  data  window  was  used.  By  using  a  tapered  rectangular  window, 

though,  I  obtained  the  reconstructed  object  of  Figure  29b.  I  noticed  that  if  a  certain 

(N  N  N  \ 
7,7,  j,  etc.j, 

then  it  would  not  reconstruct  at  all,  either.  Figure  30a  is  an  example  of  this  object 
and  Figure  30b  is  its  reconstructed  shape  after  windowing  by  a  Bartlett  window.  In 
all  cases  where  an  indeterminate  reconstruction  arose,  data  windowing  allowed  for 
its  reconstruction  (though  the  window  function  is  visible  in  the  reconstruction). 

Finally,  in  the  1-D  case,  objects  placed  in  very  large  data  records  (1X256,  512, 
1024,  etc.)  would  not  reconstruct  but  give  an  indeterminate  condition  unless  data 
windows  were  used. 

I  attempted  to  use  data  windowing  in  the  2-D  reconstruction  but  was  never  suc¬ 
cessful  in  removing  the  sine  blur.  Thus,  I  conclude  that  the  cause  of  the  sine  blur  is 
more  fundamental  in  the  2-D  case  then  mere  aperiodic  sampling  errors.  It  appears 
that,  since  my  2-D  phase  reconstruction  algorithm  only  samples  the  two  lowest 
spatial  frequency  bispectral  subplanes  [(1,0), (0,1)],  I  have  bandlimited  the  result 
to  these  lower  spatial  frequencies.  Higher  spatial  frequencies  present  in  the  recon- 
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strutted  image  would  give  me  sharper  edges  and  not  the  sine-blurred  edges  I  am 
getting  now. 


Figure  28.  Bispectrallv  Reconstructed  Padded  Tribar  with  Data  Window 
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Figure  29b.  Bispectral  Reconstruc¬ 
tion  of  (a)  Using  a  Data  Window 


Figure  30b.  Bispectral  Reconstruc¬ 
tion  of  (a)  Using  a  Data  Window 


5.2  Speckled 

To  begin  the  speckled  section  I  should  review  at  the  outset  the  flow  of  the  process¬ 
ing  that  I  performed  to  reconstruct  a  speckled  object.  Before  I  do  that  I  will  briefly 
discuss  it.  The  overall  methodology  will  be,  in  general,  the  same  whether  using 
thel-D  or  2-D  bispectrum  algorithm. 

The  truth  object  is  created  and  stored  in  an  array  that  does  not  change.  For  a 
selected  number  of  images,  M,  each  one  is  speckled  independently  using  one  of  the 
speckle  algorithms.  The  speckled  object  varies  from  frame-to-frame  by  using  a 
different  random  number  seed  to  create  the  speckle  for  each  frame.  I  then  ran¬ 
domly  jitter  the  speckled  object  in  x  and  y  for  each  frame.  This  models  image 
capture.  I  then  perform  a  homomorphic  transformation  before  I  Fourier  transform 
the  image.  I  then  calculate  the  bispectrum  of  this  function  and  sum  it  with  previ¬ 
ously  calculated  values.  Once  all  images  have  been  summed,  I  normalize  the  resul¬ 
tant  function  and  calculate  the  Fourier  magnitude  and  phase  from  their  respective 
1-D  or  2-D  algorithms.  Combining  these  quantities  and  inverse  transforming,  then 
exponentiating  that  result,  gives  the  final  answer.  A  median  filter  is  included,  if 
desired.  See  Appendix  C  for  the  flow  diagram  of  this  entire  procedure. 

As  noted  previously  several  times,  the  MLE  of  an  undegraded  version  of  the  speck¬ 
led  image  is  the  ensemble  average  of  the  multiple  speckled  images.  To  find  the 
MLE  of  the  undegraded  image  using  multiple  observations  of  the  multiplicative 
noise  model,  1  can  write: 

gk(m,n)  =  I(m,n)uk(m,n)  (81) 
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where  the  values  take  the  same  meaning  as  Equation  34  except  that  it  is  the  k 
estimate  for  N  observations  (k=l,2,3 . N).  Performing  the  homomorphic  trans¬ 

formation,  I  get: 

ln(gk(m,n))=  ln(I(m,n))  +  ln(uk(m,n))  (82) 

From  Raghuveer,  let  these  values  equal,  respectively,: 

Zk  =  6  +  nk  (S3) 

To  find  the  MLE,  I  want  to  maximize  the  probability  that  the  N  observations  pro¬ 
duce  the  original  object.  For  N  independent  observations  I  can  write: 

P(z|e)  =p(zi|e)p(z2|e  )p(z3|e  ) . p(zje)  (S4) 

where  z  is  a  vector.  Looking  only  at  one  observation,  I  can  express  it  in  terms  of 
nj  =  zx  -  0  as: 

P(zj0)  =  P(nj)  =p(z1-  0  )  (85) 


From  the  pdf  for  the  logarithm  of  the  intensity,  1  can  be  write: 


P  (  z j  —  0  )  =  exp  [-(z^  6  )  -exp[-(zj-  0)]] 


(86) 


therefore,  for  all  N  observations,  taking  the  logarithm  of  Equation  86  gives: 


In  (p(  zk-  0  ))  .  E(-(zk-8)-exp[-(zk-  9)]  )  (S') 
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(88) 


ln(p(  z.-  6 ))  =  Iz.  +  N6  -  Z  exp  [-(z.-0)] 

K  k=l  K  k=l  K 

differentiating  this  with  respect  to  0  and  setting  it  equal  to  zero  gives: 

N  =  exp  [0]  Z  exp  [  -  z  ]  (89) 

k=l  K 

returning  to  the  original  definitions,  I  can  substitute  in  for  0  and  zk  and  divide  by  N 
to  get: 


I(  m,  n  )  Z  gk ( m,  n  )  (90) 

N  k=i  k 

Thus,  averaging  N  observations  of  the  noisy  image  will  provide  the  MLE  of  the 
degraded  image.  This  result  is  similar  to  the  result  given  by  Lim  and  Nawab  (1981) 
except  they  kept  the  noise  expression  in  the  non-logarithmic  domain  during  the 
proof. 

To  see  what  this  implies  for  image  reconstruction,  notice  the  speckled  version  of  a 
sine2  seen  in  Figures  31a  and  b,  both  in  2-D  and  1-D,  respectively.  Figures  32a 
and  b  are  the  resultant  images  after  150,  independently-speckled  images  were 
ensemble  averaged.  Notice  that  you  do  not  reconstruct  the  image  to  obtain  Figure  9 
but,  hopefully,  to  see  what  the  object  was  in  the  first  place!  This  is  a  result  of  the 
severe  image  degradation  caused  by  the  fully-developed  speckle. 


-98- 


31b.  One-Dimensional  Slice  of  (a)  DDM  Used 
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ure  32h.  One-Dimensional  Slice  of  (a) 
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5.2.1  One-Dimensional  Examples 


The  first  example  used  showing  image  reconstruction  is  found  in  Figures  33a-g. 
Figure  33a  is  the  truth  object  and  Figure  33b  is  a  typical  speckled  object  before  the 
bispectrum  routine  begins.  The  DDM  was  used  to  speckle  the  object.  Figure  33c  is 
the  ensemble  averaged  reconstruction  of  the  jittered  images  and  Figure  33d  is  the 
reconstructed  object  after  200  images  were  bispectrum  averaged.  Figure  33e  is  the 
object  after  a  1X3  median  filter  has  passed  over  the  reconstructed  image.  Figure 
33f  is  the  reconstructed  object  if  no  homomorphic  transformation  is  taken.  As  one 
can  see,  no  object  reconstruction  will  occur  for  speckled  objects  unless  the  natural 
logarithm  is  taken  first.  Figure  33g  is  the  reconstructed  object  after  only  30  images 
were  bispectrum  averaged,  showing  how  rapidly  bispectral  averaging  can  converge 
to  a  solution. 

One  reason  the  bispectrum  rapidly  converges  to  a  solution  is  shown  in  Figures 
34a-c.  Figure  34a  is  the  natural  logarithm  of  a  sinc2(x/8)  function.  Wherever  the 
real  function  goes  to  zero  I  set  its  density  domain  representation  to  -150.  Figure 
34b  is  the  log  of  the  speckle  intensity  shown  in  Figure  33b.  Notice  how  the  noise 
now  appears  to  “ride"  on  top  of  the  function,  essentially  improving  the  SNR  of  the 
resultant  image.  Figure  34c  is  the  bispectrally  averaged  estimate  of  the  speckled 
image  prior  to  exponentiation  (which  then  gives  Figure  33d).  Notice  how  the  noise 
has  smoothed  considerably. 
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Figure  33a-d.  Example  of  Bispectral  Averaging  of  Jittered,  Speckled  Objects 
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(g) 

Figure  33e-g.  Example  of  Bispectral  Averaging  of  Jittered,  Speckled  Objects 

(Corn'd) 
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Figure  34a-c.  Result  of  Bispectral  Averaging  an  Object  in  its  Logarithmic 

Domain 


Another  reason  the  bispectrum  can  handle  speckled  images  is  seen  in  Figures 
35a-f.  I  use  the  same  sine2  again  but  with  no  zero  padding  so  as  to  view  the  results 
better.  I  show  the  logarithmic  representation  but  I  also  include  the  Fourier  trans¬ 
form  of  the  log  of  the  speckle  intensity  (Figure  35f).  Compare  this  to  Figure  35e. 
Thus,  the  bispectrum,  which  uses  the  function  in  Figure  35f  to  obtain  its  values, 
basically  sees  a  very  close  approximation  of  the  original  object  when  processing  in 
this  domain  (the  Fourier  of  the  log  domain). 

The  next  three  sets  of  Figures  use  a  tri(x)  object  to  describe  three  key  points. 
Figures  36a-f  show  a  non-jittered  object  with  its  bispectrally-averaged  rendition. 
Notice  that,  even  using  the  ensemble  average  (Figure  36c)  I  only  get  a  form  of  the 
tri  back.  Figure  36d  is  what  a  1X3  median  filter  would  do  if  used  directly  on 
Figure  36b.  Figures  36e-f  are  the  bispectrally  reconstructed  image  without  and 
with  the  median  filter,  respectively. 

The  next  set  of  Figures  (37a-e)  are  of  a  tri  that  is  jittered  randomly  in  a  wider  field 
200  times.  Since  averaging  is  done  in  the  bispectral  domain,  the  shifting  is  ignored 
by  the  bispectrum  and  reconstruction  occurs  (Figures  37c-d).  Figure  37e  shows  the 
result  if  no  homomorphic  transformation  occurs. 
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Figure  35a-d.  Comparison  of  the  Fourier  Representation  of  a  Speckled  Object 

with  its  Unspeckled  Version 
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Figure  35e-f.  Comparison  of  the  Fourier  Representation  of  a  Speckled  Obje 

with  its  Unspeckled  Version  (Corn'd) 
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Figure  36e-f.  Example  of  Bispectral  Averaging  of  Speckled  Object  (Corn'd) 
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Figure  37a-c.  Example  of  Bispectral  Averaging  of  Jittered,  Speckled  Object 
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Figure  37d-e.  Example  of  Bispectral  Averaging  of  Jittered,  Speckled  Object 

(Cont’d) 
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The  third  set  of  Figures  (38a-e)  show  why  bispectral  averaging  works  for  the  tri. 
The  log  of  the  speckled  tri  (Figure  38b)  has  such  an  improved  SNR  (visually)  than 
the  speckled  version.  After  averaging  (Figure  38c),  the  noise  is  considerably 
smoothed.  Figure  38d  is  the  Fourier  transform  of  the  log  of  the  tri.  Compare  this 
to  Figure  38e:  there  is  no  difference!  For  the  scale  factors  used,  the  noise  is  practi¬ 
cally  nonexistent  in  this  domain.  This  is  what  the  bispectrum  sees  when  it  is  being 
calculated. 

Finally,  Figures  39  and  40  are  speckled  and  reconstructed  versions  of  Figures  18a 
and  22a,  respectively.  Notice  how  well  we  can  discriminate  the  form  of  these  ob¬ 
jects  which  was  lost  in  the  speckle  (see  Korwar  and  Pierce,  (1981)  for  a  discussion 
of  this  loss  of  discrimination)  after  bispectral  averaging. 

5.2.2  Two-Dimensional  Examples 

I  was  so  pleased  with  the  results  of  the  1-D  algorithm  that  I  decided  to  extend  it  to 
2-D  images.  To  do  this,  I  would  raster-scan  the  2-D  image,  after  jittering  and 
speckling,  to  reconstruct  the  object.  Before  attempting  this,  though,  I  tested  the 
1-D  algorithm  on  Figure  41a,  which  is  the  object  “E”  jittered  in  a  “smearing” 
action  (Figure  41b).  Figure  41c  is  the  result:  perfect  reconstruction! 

I  then  tried  it  on  a  very  small  amount  of  additive  Gaussian  noise.  It  is  so  small 
(variance  +  0.003  about  the  image)  that  you  can  barely  see  it  on  a  3-D  plot  (Figure 
42a).  After  32  iterations,  an  “E”  results  that  is  slightly  distorted  (Figure  42b).  By 
150  iterations,  the  noise  is  gone  and  the  original  image  results  (though  it  is  phase- 
shifted)  (Figure  42c). 
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3Sa-c.  Bispectral  Reconstruction  in  Logarithmic  Domain  and  Fourier 
Comparison  of  Speckled  Object  with  Original  Version 
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Figure  39a-b.  Bispectral  Reconstruction  of  Speckled  Object 
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Figure  41a-c.  Bispectral  Reconstruction  of  2-D  Object  Using  1-D 

Algorithm  -  No  Noise 
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I  next  increased  the  variance  of  the  zero  mean,  additive  Gaussian  noise  to  10X  the 
previous  value  to  see  the  results.  Figure  43b  is  the  result:  unrecognizable  recon¬ 
struction!  When  noise  was  present,  a  2-D  extension  of  the  1-D  algorithm  was  not 
possible.  This  is  due  to  the  phase  errors  propagating  throughout  the  reconstruction 
algorithm  from  the  recursion  formula  used.  Larger  noise  means  larger  phase 
errors. 

Obviously,  speckle  was  an  even  worse  problem  than  this  so  I  had  to  find  a  2-D 
algorithm  to  work  on  2-D  images.  That  2-D  algorithm,  developed  previously  in 
this  thesis,  was  used  on  these  2-D  images. 

Figure  44a  is  the  speckled  version  of  Figure  27a  using  the  CWM  with  a  contrast  of 
1.053862.  Figure  44b  shows  30  speckled  images  ensemble-averaged  together.  Fig¬ 
ure  44c  is  the  reconstructed  object  after  30  images  were  bispectrally  averaged,  and 
Figure  44d  is  the  same  image  but  with  a  median  filter  passed  over  the  image. 


i 
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Figure  44a-d.  Bispectral  Reconstruction  Using  2-D  Algorithm.  Object 

Speckled  Via  CWM 


Figure  45a-c.  Bispectral  Reconstruction  Using  2-D  Algorithm.  DDM  Used  to 

Speckle  Object 
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Figure  46a-d.  Bispectra]  Reconstruction  Using  2-D  Algorithm.  AGM  Used 

Speckle  Object 


Figure  47  is  the  reconstructed  spacecraft  from  30  speckled  images  using  a  correla¬ 
tion  technique  I  devised  to  keep  track  of  the  shift  of  the  object  from  frame-to- 
frame.  Basically,  I  use  the  first  speckled  image  as  a  reference  with  which  I  corre¬ 
late  all  of  the  following  images.  Since  speckle  is  wideband,  random  noise,  their 
correlations  have  a  distinct  peak  when  they  overlap  during  the  correlation  proce¬ 
dure.  I  keep  track  of  this  peak’s  location  and  align  the  images  with  the  reference 
image  so  that  ensemble  averaging  can  be  performed.  The  2-D  contour  plot  shows 
a  gradual  sloping  toward  the  center  of  each  object  structure  and  the  edges  are  a  bit 
blurred.  No  homomorphic  transformation  is  taken  with  this  method. 

A  very  important  point  has  to  be  made  here.  When  using  the  bispectrum  to  esti¬ 
mate  a  speckled  object,  the  logarithm  must  be  taken  prior  to  the  estimation.  This 
was  seen  in  the  1-D  case  and  Figure  48  is  the  proof  in  the  2-D  case  for  a  speckled 
binary  spacecraft. 

In  using  the  homomorphic  transformation,  I  also  discovered  that  care  must  be 
taken  prior  to  exponentiation  to  obtain  optimum  results.  In  some  cases,  the  recon¬ 
structed  Fourier  object  had  values  that  would  cause  an  arithmetic  overflow  if  not 
normalized  by  some  constant.  The  choice  of  this  constant  also  was  an  integral  part 
of  the  quality  of  the  reconstruction.  Different  objects  required  different  normaliz¬ 
ing  factors  and  the  choice  was  purely  subjective  based  on  what  the  reconstructed 
object  needed  to  enhance  its  appearance,  i.e,  edges.  Since  the  goal  of  speckle-im¬ 
age  reconstruction  is  the  return  of  the  object  form  destroyed  by  the  speckle,  this  is 
the  only  necessary  metric  by  which  to  base  the  decision  on  the  choice  of  the  nor¬ 
malizing  factor. 
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5.3  Conclusions 


This  thesis  has  shown  that,  conceptually,  bispectrum  estimation  is  robust  enough  to 
recover  speckle-degraded  images  when  a  homomorphic  transformation  is  taken 
and  the  objects  are  averaged  in  the  bispectrum.  This  is  especially  true  in  the  1-D 
case,  where  any  object  could  be  reconstructed  from  any  form  of  the  speckle  model. 
The  necessary  object  conditions  were  that:  1)  it  had  to  stay  consistently  within  the 
field  of  view  of  the  array  when  jittering  and,  2)  even  though  it  could  change  posi¬ 
tion  within  the  field  of  view,  it  could  not  change  orientation. 

The  phase  recovery  problems  with  the  2-D  case  were  not  totally  resolved.  Phase 
unwrapping  was  thought  to  be  a  solution  but  trying  to  unwrap  a  2-D  object’s  phase 
to  obtain  a  unique  solution  in  2-D  is  difficult,  especially  if  the  phase  changes  by 
modulo  27r  every  pixel!  A  2-D  object  has  phase  terms  going  in  the  x  and  y  direc¬ 
tions;  phase  unwrapping  in  one  direction  may  not  be  what  is  needed  in  the  other 
direction  at  each  m,n  pixel  location.  This  is  because  phase  unwrapping  attempts  to 
impose  continuity  on  the  phase  function  to  remove  any  ambiguities  present  in  the 
phase  (Fiddy  and  Stark).  One  ambiguity  that  arises  is  from  the  introduced  phase 
shift  of  the  reconstructed  bispectrum.  If  I  could  write  the  reconstructed  estimate  as 
an  integer  phase  shift  then  I  would  get  (after  Kaveh  and  Soumekh): 

F  (  k  £  )  =  F(  X  5  )  exp  [  j  2 77  n  ]  =  |f  (  k  £  )  |exp  [j  ( <J>p(  X  t  )  +  2ir  n)]  (91) 
The  argument  of  F  ( k,  £  )  for  the  phase,  then,  is  given  by: 

ARG{F  (  k,  £) }  =  27in  +  arctan  (  —  ■  ^ )  1  (92) 

L  J  VKe(F(\,£))J 
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The  ambiguity  for  the  value  of  the  integer,  n,  is  no  problem  since 
<1>F ( X •  ^ )  and  $F(\,£  )  +  2Trn  give  the  same  function  for  the  given  modulus, 

j  F(X.£  )  |.  The  term  2im  carries  no  information  thus  no  physical  means  exist  to 
determine  it  anyway  (Takajo  and  Takahashi,  1988).  But  if  I  take  the  logarithm  of 
the  function,  F(X,£),  we  need  to  know  n.  This  can  be  seen  from  the  logarithm  of 

the  complex  function,  F  ( M  £ ),  which  is  written  as: 

In  (F  (  X,  £  )  )  =  In  (j  F(  X,  £  )  |  }  +  j  +  arctan(^(^(  ^  [  ) )))  (93) 

Note  that  the  log  of  the  reconstructed  object  has  the  same  real  part  as  the  log  of 
the  real  object  but  it  has  an  imaginary  part  that  differs  by  integrals  of  2t r.  The 
logarithm  of  the  spectrum  is  no  longer  unique,  as  in  the'  real  object’s  case,  due  to 
this  phase  ambiguity.  Phase  unwrapping  becomes  a  process  whereby  we  attempt  to 
retrieve  this  term. 

Finally,  phase  can  be  described  physically  as  the  wavefront  surface  of  a  propagat¬ 
ing  complex  wave  field.  If  the  amplitude  of  the  propagating  wave  is  zero  at  any 
point,  its  phase  is  indeterminate.  Moving  about  that  point  causes  a  phase  change  of 
2tt  (1  wavelength)  implying  wavefront  dislocations  (Fiddy  and  Stark).  Since 
speckle  is  caused  by  modulo  2tt  phase  excursions,  there  are  as  many  wavefront 
dislocations  as  there  are  speckles.  Phase  unwrapping  in  this  case  would  be  indeter¬ 
minate  since  the  chance  of  obtaining  an  unambiguous  phase  unwrap  would  be 
slim. 
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In  summation,  this  thesis  was  an  attempt  to  apply  spectral  techniques  to  remove 
multiplicative  noise  from  an  object  jittered  against  a  uniform  background  in  both 
one  and  two  dimensions.  It  was  also  an  exercise  in  determining  the  utility  of  the 
bispectrum  as  a  signal  processing  tool.  The  homomorphic  transformation  creates  a 
situation  where  you  have  a  signal  plus  noise  and  this  noise  is  small  compared  to 
the  signal  in  the  bispectral  domain.  Averaging  in  this  domain  allows  us  to  obtain 
an  estimate  of  the  original  object.  Finally,  this  thesis  attempts  to  review  speckle 
theory  and  then  model  it  in  different  fashions  to  fully  understand  the  mechanisms 
behind  coherent  speckle.  Further  work  could  attempt  to  use  actual  speckled  images 
and  also  increase  the  robustness  of  the  2-D  phase  reconstruction  routine.  One 
could  also  look  at  using  additional  enhancement  routines  (aside  from  the  median 
filter  used  in  this  thesis)  as  post-processing  techniques  to  improve  image  recon¬ 
struction. 
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6.0  APPENDICES 


Appendix  A 

Tabular  Synopsis  Of  Speckle  Models 


Model 

Description 

MMMMi-t*  1 1  n  tj « m  i  tJBMSi 

Transfer  Function? 

DDM 

Circular,  complex 

Gaussian  r.v.  with  zero 
mean,  i.i.d.  real  and 
imaginary  components 
created  as  field. 

Magnitude  squared  for 
speckle  intensity. 

Large  number  of 
surface  scatterers; 
no  multiple  scat¬ 
tering;  no  depol¬ 
arization  upon 
reflection. 

NO-assumed  embedded 
in  incoherent  image. 

TFM 

Speckle  field  as  in  DDM, 
then  Fourier  transformed 
for  propagation.  Coherent 
transfer  function  multiplied 
by  this  and  then  inverse 
transformed  before 
magnitude  squaring. 

Coherent  optical 
system  can  resolve 
all  details  in  object. 

YES-user  defined.  A 
key  player  in  ultimate 
object  appearance  and 
statistics. 

CWM 

Speckle  field  as  in  DDM. 
Complex  spatial  window 
of  equal  strengths 
correlates  complex  field. 
Magnitude  squaring 
result  gives  intensity. 

High  quality  optics 
w'ould  correlate 
the  speckle  locally. 

YES-in  the  form  of  a 

3X3,  5X5,  etc.  spatial 
window.  Window  acts 
as  averaging  filter. 

RPM 

Modulo  2tr  random  phase 
created  from  uniform 
random  number  generator. 
This  used  with  image 
field  via  Euler’s  notation 
to  create  speckle  field, 
then  magnitude  squared 
to  give  speckle  intensity. 

Can  be  imaged  or 
pupil-plane  speckle. 

Object  surface  is 
rough  on  the  order 
of  the  wavelength 
of  the  incident 
coherent  source, 
up  to  multiples  of 

2tt  until  coherence 
length  of  source  is 
reached. 

YES-is  optional  but  would 
be  implemented  similar  to 

TFM  once  speckle  field  has 
been  created. 

NEM 

Noise  function  created 
that  has  a  negative 
exponential  pdf.  This  is 
then  multiplied  by  the 
incoherent  image  of 
the  truth  object. 

Speckle  noise  is  a 
signal-dependent 
and  multiplicative 
noise  process. 

NO-field  statistics  are  never 
derived. 

AGM 

Zero  mean,  unity  variance, 
Gaussian  noise  added  to 
the  logarithm  of  the 
object  intensity. 

Log  of  multiple 
speckle  intensity 
is  approximately 
Gaussian. 

NO-field  statistics  are  never 
derived. 
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Appendix  B 

Bispectrum  Fourier  Expression 

This  appendix  is  proof  of  Equation  43  in  the  text.  Rewriting  Equation  42,  the 
bispectrum  is  given  as  the  Fourier  transform  of  the  triple  correlation: 

oo  oo 

I(3)(fj,  f2)  -  f  Ji(3Wx2)exp[-27rj(f1x1+f2x2)]  dxj  dx2  (Bl) 
-00-00 

oo 

where:  i(3)(xrx2)=Ji(x)i(x  +  xt)i(x  +  x2 )  dx  (B2) 

-OO 

is  the  triple  correlation  of  the  spatial  image  i(x)  and  x  is  a  2-D  vector.  Substituting 
B2  into  Bl  gives  the  bispectrum  definition  in  terms  of  the  real  object,  i(x)  and  its 
shifted  versions,  i(x  +  Xj )  and  i(x  +  x2)  as: 

oo  oo  oo 

I<3>  ( f  i ,  f2)  -/  J  /i(x)i(x +Xj  )i(x +x2) 

.00  —  00  — oo  (B ') 

exp  [-2iTj(f1  x  j  +  f2  x 2  )  ]  tlx  dxj  dx2 


let  (3  =  (x  +  x1),a  =  (x+  x2),  therefore  d3  =  dxj  and  da  =  dx2  since  x  is 
constant  and  (Xj+  x2)are  the  shifted,  changing  values.  Substitute  these  into  B3 
above  to  get: 


oo  OO  OO 

-  /  /  /i(x)i(p)i(a) 

—00  —  00  -oo 

exp  [-2tt  j(fj  ( 3  -  x )  +  f2  (a  -  x  ) )]  dx  d3da 


(B4) 
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00 


Collecting  terms,  noting  that  J  f(x)  exp  [-  j2Trj;x]  =F(£)  is  the  continuous 

-OO 

definition  of  the  Fourier  transform  of  f(x)  (Gaskill),  I  get: 


9? 


I(3)(f  f  )  =  /i(x)dxJiO)  exp[-j  277f13]d3  Ji(a) 

—  OO  —  OO  —  oo 

exp  [-j  2tt  f2  a]  da  exp  [  j  2tt fj  x]  exp  [  j  2'rr  f2x] 


(B5) 


l(3)(fj.  f2)  =fi(*)dx  exp  [j2irx(f1  +  f2)]  I(f  )I(f ,)  (B6) 

—  OO 

Expanding  the  exponential  term  in  Equation  B6  into  its  Euler’s  notation,  exp[-j0]  = 
cos(0)  -  jsin(0),  I  can  write  the  Fourier  transform  definition  in  terms  of  its  real  and 

imaginary  components,  respectively,  as  F(0  *  FRe(£)  ~  jFim(0.  Thus,  I  can 
write  Equation  B6  as: 

l‘3)(fl.  f2)-I(fl)I(f2)[lRe(fl+  f2)  +  j  Ilm(fl+  f2)] 

OO 

Noting  that  F  *  (£)  =  Jf(x)  exp  [ j 2'T7'^x  ]  Equation  B7  then  becomes: 

-OO 

I(3>(fi.  f2) -l(fi)I(f2)I  *Ui  +  f2) 


(B7) 


(BS) 
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From  this  expression  I  can  simply  write  the  bispectrum  magnitude  and  phase  rela¬ 
tionships  as  functions  of  the  object’s  Fourier  representation  as  (  |  I  |  is  the  object 
magnitude  and  4>  is  the  object  phase): 


I(3)(fl,  f 2) 

=  |l(3)(fi.  h)  |  exp  [  jt|r  (flt  f2)] 

09) 

where:  |l(3j  (fi,  f2) | 

=  |  I(fl)|  |  I  (f  2 )  |  |l(-fl-f2)| 

(BlOa) 

*  (ft.  f2) 

=  4>(fi)  +  4>(f2)  -  +  f2) 

(BlOb) 

To  compare  this  result  with  the  power  spectrum  density  (PSD), 

Fourier  transform  of  the  autocorrelation,  it  is  defined  as: 

which  is  the 

i<2)(fi)  = 

l  i(2>(fi)i  «p  [ j ^  (fi)] 

(Bll) 

where:  1 1 ( 2 '  (  fi )  1  = 

h(fi)l  iK-f,)i  - 1  Kfi)p 

(B12a) 

\|/(fl) 

4>(fi)  -<l>(fi)  -  o 

(B12b) 

The  PSD  gives  us  magnitude  information  regarding  the  object  but  phase  is  lost. 
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Appendix  C 

Bispectrum  Reconstruction  Algorithm  Flowchart 

This  top  level  block  diagram  gives  the  logic  flow  of  how  the  algorithm  was  imple¬ 
mented  to  perform  shift-invariant,  speckle-degraded  image  reconstruction.  The 
details  of  each  block  are  in  the  text. 
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(RANDOMLY  JITTER  THE  SPECKLED  OBJECT! 

li  — ) IMAGE  CAPTURE  UOOELED 
[TAKE  NATURAL  LOG:  h(m,n)=ln(g(m,n))l 

ii 

(TAKE  FOURIER  TRANSFORM:  H(i,l)=3{h(m,n))| 
[CALCULATE  AND  SUM  BISPECTRUM  SUBPLANESl 

li 

_ YES  ANOTHER 

IMAGE? 

li  NO 

[AVERAGE  THE  COMPUTED  BISPECTRUMl 

li 

[CALCULATE  BISPECTRUM  MAGNITUDE  AND  PHASE! 

, _ 4 . . 

ICOMPUTE  FOURIER  MAGNITUDE  FROM  I  B  l[ 

li 

ICOMPUTE  FOURIER  PHASE  FROM  ^1 

li 

[COMBINE  COMPONENTS  INTO  COMPLEX  FUNCTION! 

li 

IlNVERSE  TRANSFORM  COMPLEX  IMAGE! 

li  -►  DIVISOR? 

[EXPONENTIATE  INVERSE-TRANSFORMED  IMAGE] 

li 

MEDIAN  FILTER? 

li 

[QUANTIZE  AND  DISPLAY  RECONSTRUCTED  IMAGE! 


Figure  C-l.  2-D  Image  Reconstruction  Overview 
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